In Equation 3.6, $\rho$ is a block diagonal matrix and $u$ is a realization of a random vector $\mathcal{U}$. The diagonal elements of $\rho$, say $\rho$, are matrices which contain the correlations between $\gamma_1$ and $\gamma_2$. For example, if $\gamma_1 = (\gamma_1^I, \gamma_1^S)^\top$ and $\gamma_2 = (\gamma_2^I, \gamma_2^S)^\top$, with $I$ = Intercept and $S$ = Slope,

$$ \rho = \begin{pmatrix} \operatorname{corr}(\gamma_1^I, \gamma_2^I) & \operatorname{corr}(\gamma_1^I, \gamma_2^S) \\ \operatorname{corr}(\gamma_1^S, \gamma_2^I) & \operatorname{corr}(\gamma_1^S, \gamma_2^S) \end{pmatrix} \quad \text{and} \quad \rho = \operatorname{diag}(\rho, \dots, \rho). \qquad (3.7) $$

The bivariate linear mixed-effects model is then re-expressed as:

$$ \begin{aligned} y_1 &= X_1\beta_1 + Z_1\Lambda_{\theta_1}u_1 + \varepsilon_1, \\ y_2 &= X_2\beta_2 + Z_2\Lambda_{\theta_2}u_2 + \varepsilon_2, \end{aligned} \qquad (3.8) $$

with

$$ u = \begin{pmatrix} u_1 \\ u_2 \end{pmatrix} \sim \mathcal{N}(0, \Sigma_u), \quad \varepsilon = \begin{pmatrix} \varepsilon_1 \\ \varepsilon_2 \end{pmatrix} \sim \begin{pmatrix} 0 & \sigma_1^2 \mathbf{I}_N \\ 0 & \sigma_2^2 \mathbf{I}_N \end{pmatrix}. \qquad (3.9) $$

Then the parameters which will be estimated are $\beta_1, \beta_2, \sigma_1^2, \sigma_2^2, \theta_1, \theta_2$ and $\rho$.

## 3.3 Parameters' estimates

In this Section, we first provide the likelihood of the model's parameters and then give the REML criterion which will be optimized for the obtaining of the parameters' REML estimates.

### 3.3.1 ML criterion

The ML criterion is the log-likelihood of the model's parameters which is displayed through the following theorem

**Theorem 3.3.1.** Suppose that $\mathbf{y} = (y_1^\top, y_2^\top)^\top$ satisfies the bivariate linear mixed-effects model expressed by Equations (3.8 and 3.9), where $\beta_1, \beta_2, \sigma_1^2, \sigma_2^2, \theta_1, \theta_2, \rho$ are the parameters which need to be estimated, and $\beta = (\beta_1^\top, \beta_2^\top)^\top$, $\sigma = (\sigma_1^2, \sigma_2^2)^\top$, $\theta = (\theta_1^\top, \theta_2^\top)^\top$. Denoting by $Y_\sigma = (\sqrt{\sigma_2^2}y_1^\top, \sqrt{\sigma_1^2}y_2^\top)^\top$, $X_\sigma = (\sqrt{\sigma_2^2}X_1 \quad 0 \quad 0 \quad \sqrt{\sigma_1^2}X_2)$, $Z_{\sigma\theta} = (\sqrt{\sigma_2^2}Z_1\Lambda_{\theta_1} \quad 0 \quad 0 \quad \sqrt{\sigma_1^2}Z_2\Lambda_{\theta_2})$, and $\mu_{\mathcal{U}|\mathbf{y}=\mathbf{y}$ the conditional mean of $\mathcal{U}$ given that $\mathcal{Y} = \mathbf{y}$, the log-likelihood of $\beta, \sigma, \theta$ and $\rho$ given $\mathbf{y}$ is expressed as

$$ \ell(\beta, \theta, \rho, \sigma | \mathbf{y}) = - \frac{r(\hat{\beta}_{\theta, \rho, \sigma}, \mu_{\mathcal{U}|\mathbf{y}=\mathbf{y}}) + \|R_X(\beta - \hat{\beta}_{\theta, \rho, \sigma})\|^2}{2\sigma_1^2\sigma_2^2} - \frac{N-q}{2} \log(\sigma_1^2\sigma_2^2) \\ - \frac{1}{2} \log(|\Sigma_u|) - \frac{1}{2} \log(|L_{\theta, \rho, \sigma}|^2), \qquad (3.10) $$
where $q = q_1 + q_2$, $\hat{\beta}_{\theta,\rho,\sigma}$ and $\mu_{\mathcal{U}|\mathcal{Y}=\mathbf{y}}$ satisfy

$$ \begin{pmatrix} X_{\sigma}^{\top} X_{\sigma} & X_{\sigma}^{\top} Z_{\sigma\theta} \\ Z_{\sigma\theta}^{\top} X_{\sigma} & Z_{\sigma\theta}^{\top} Z_{\sigma\theta} + \sqrt{\sigma_1^2 \sigma_2^2 \Sigma_u^{-1}} \end{pmatrix} \begin{pmatrix} \hat{\beta}_{\theta,\rho,\sigma} \\ \mu_{\mathcal{U}|\mathcal{Y}=\mathbf{y}} \end{pmatrix} = \begin{pmatrix} X_{\sigma}^{\top} \\ Z_{\sigma\theta}^{\top} \end{pmatrix} Y_{\sigma}, \quad (3.11) $$

$$ r(\hat{\beta}_{\theta,\rho,\sigma}, \mu_{\mathcal{U}|\mathcal{Y}=\mathbf{y}}) = \left\| Y_{\sigma} - X_{\sigma} \hat{\beta}_{\theta,\rho,\sigma} - Z_{\sigma\theta} \mu_{\mathcal{U}|\mathcal{Y}=\mathbf{y}} \right\|^2 + \sigma_1^2 \sigma_2^2 \mu_{\mathcal{U}|\mathcal{Y}=\mathbf{y}}^{\top} \Sigma_u^{-1} \mu_{\mathcal{U}|\mathcal{Y}=\mathbf{y}}, \quad (3.12) $$

$L_{\theta,\rho,\sigma}$ satisfies

$$ L_{\theta,\rho,\sigma} L_{\theta,\rho,\sigma}^{\top} = Z_{\sigma\theta}^{\top} Z_{\sigma\theta} + \sqrt{\sigma_1^2 \sigma_2^2 \Sigma_u^{-1}}, \quad (3.13) $$

and $R_X$ satisfies

$$ \begin{pmatrix} X_{\sigma}^{\top} X_{\sigma} & X_{\sigma}^{\top} Z_{\sigma\theta} \\ Z_{\sigma\theta}^{\top} X_{\sigma} & L_{\theta,\rho,\sigma} L_{\theta,\rho,\sigma}^{\top} \end{pmatrix} = \begin{pmatrix} R_X & \mathbf{0} \\ R_{ZX} & L_{\theta,\rho,\sigma}^{\top} \end{pmatrix}^{\top} \begin{pmatrix} R_X & \mathbf{0} \\ R_{ZX} & L_{\theta,\rho,\sigma}^{\top} \end{pmatrix}. \quad (3.14) $$

*Proof.* Denoting by $f_X(\cdot)$ the density function of any random vector $\mathcal{X}$,

$$ f_{\mathcal{Y}}(\mathbf{y}) = \int_{\mathbb{R}^{q_1+q_2}} f_{\mathcal{Y}}(\mathbf{y}, \mathbf{u}) d\mathbf{u}, \quad (3.15) $$

where

$$ \begin{aligned} f_{\mathcal{Y}}(\mathbf{y}, \mathbf{u}) &= f_{\mathcal{Y}|\mathcal{U}}(\mathbf{y}|\mathbf{u}) f_{\mathcal{U}}(\mathbf{u}) = f_{Y_1|U_1}(y_1|u_1) f_{Y_2|U_2}(y_2|u_2) f_{\mathcal{U}}(\mathbf{u}) \\ &= (2\pi\sigma_1^2)^{-N/2} (2\pi\sigma_2^2)^{-N/2} (2\pi)^{-q_1+q_2/2} |\Sigma_{\mathbf{u}}|^{-1/2} \exp\left(-\frac{\|y_1 - X_1\beta_1 - Z_1\Lambda_{\theta_1}u_1\|^2}{2\sigma_1^2}\right. \\ & \qquad \left. - \frac{\|y_2 - X_2\beta_2 - Z_2\Lambda_{\theta_2}u_2\|^2}{2\sigma_2^2} - \frac{1}{2}\mathbf{u}^{\top}\Sigma_{\mathbf{u}}^{-1}\mathbf{u}\right). \end{aligned} \quad (3.16) $$

Let us denote by $\tilde{\Sigma}$ the matrix such that

$$ \Sigma_u^{-1} = \tilde{\Sigma}^{\top} \tilde{\Sigma}. \quad (3.17) $$

It then comes that $\mathbf{u}^{\top}\Sigma_u^{-1}\mathbf{u} = \|\tilde{\Sigma}\mathbf{u}\|^2$ and

$$ \begin{aligned} & \frac{\|y_1 - X_1\beta_1 - Z_1\Lambda_{\theta_1}u_1\|^2}{\sigma_1^2} + \frac{\|y_2 - X_2\beta_2 - Z_2\Lambda_{\theta_2}u_2\|^2}{\sigma_2^2} + \mathbf{u}^{\top}\Sigma_u^{-1}\mathbf{u} \\ &= \frac{\|\sqrt{\sigma_2^2}(y_1 - X_1\beta_1 - Z_1\Lambda_{\theta_1}u_1)\|^2 + \|\sqrt{\sigma_1^2}(y_2 - X_2\beta_2 - Z_2\Lambda_{\theta_2}u_2)\|^2 + \|\sqrt{\sigma_1^2\sigma_2^2}\tilde{\Sigma}\mathbf{u}\|^2}{\sigma_1^2\sigma_2^2} \\ &= \left\| \begin{pmatrix} \sqrt{\sigma_2^2}y_1 \\ \sqrt{\sigma_1^2}y_2 \\ \mathbf{0}_{q_1+q_2} \end{pmatrix} - \begin{pmatrix} \sqrt{\sigma_2^2}X_1 & \mathbf{0}_{N_{p_2}} & \sqrt{\sigma_2^2}Z_1\Lambda_{\theta_1} & \mathbf{0}_{N_{q_2}} \\ \mathbf{0}_{N_{p_1}} & \sqrt{\sigma_1^2}X_2 & \mathbf{0}_{N_{q_1}} & \sqrt{\sigma_1^2}Z_2\Lambda_{\theta_2} \\ \mathbf{0}_{q_1+q_2,p_1+p_2} & & \sqrt{\sigma_1^2\sigma_2^2}\tilde{\Sigma} & \end{pmatrix} \begin{pmatrix} \beta \\ \mathbf{u} \end{pmatrix} \right\|^2 \end{aligned} \quad (3.18) $$
$$= \left\| Y_{\Lambda} - Z_{X\Lambda} \begin{pmatrix} \beta \\ u \end{pmatrix} \right\|^2 \qquad (3.19)$$

$$= g(\beta, u, \theta, \rho, \sigma). \qquad (3.20)$$

$$\begin{pmatrix} \hat{\beta}_{\theta,\rho,\sigma} \\ \mu u | \mathbf{y}=y \end{pmatrix} = \operatorname*{arg\,min}_{u,\beta} g(\beta, u, \theta, \rho, \sigma) \iff Z_{X\Lambda}^{\top} Z_{X\Lambda} \begin{pmatrix} \hat{\beta}_{\theta,\rho,\sigma} \\ \mu u | \mathbf{y}=y \end{pmatrix} = Z_{X\Lambda}^{\top} Y_{\Lambda} \text{ (normal eq.)}, \qquad (3.21)$$

with

$$Z_{X\Lambda}^{\top} Z_{X\Lambda} = \begin{pmatrix} X_{\sigma}^{\top} X_{\sigma} & X_{\sigma}^{\top} Z_{\sigma\theta} \\ Z_{\sigma\theta}^{\top} X_{\sigma} & Z_{\sigma\theta}^{\top} Z_{\sigma\theta} + \sqrt{\sigma_1^2 \sigma_2^2} \Sigma_u^{-1} \end{pmatrix} \quad \text{and} \quad Z_{X\Lambda}^{\top} Y_{\Lambda} = \begin{pmatrix} X_{\sigma}^{\top} \\ Z_{\sigma}^{\top} \end{pmatrix} Y_{\sigma}. \qquad (3.22)$$

By setting $p = p_1 + p_2$, $dim(Z_{X\Lambda}) = (2N + q) \times (p + q)$ and $S = Im(Z_{X\Lambda})$ is a subspace of $\mathbb{R}^{2N+q}$. $Y_{\Lambda} \in \mathbb{R}^{2N+q}$ and $Z_{X\Lambda} \begin{pmatrix} \hat{\beta}_{\theta,\rho,\sigma} \\ \mu u | \mathbf{y}=y \end{pmatrix}$ is the orthogonal projection of $Y_{\Lambda}$ on $S$. Then,

$$Z_{X\Lambda} u \perp \left[ Y_{\Lambda} - Z_{X\Lambda} \begin{pmatrix} \hat{\beta}_{\theta,\rho,\sigma} \\ \mu u | \mathbf{y}=y \end{pmatrix} \right], \forall u \in \mathbb{R}^{p+q}. \qquad (3.23)$$

And $g(\beta, u, \theta, \rho, \sigma)$ can then be rewritten as:

$$g(\beta, u, \theta, \rho, \sigma) = \left\| Y_{\Lambda} - Z_{X\Lambda} \begin{pmatrix} \beta \\ u \end{pmatrix} + Z_{X\Lambda} \begin{pmatrix} \hat{\beta}_{\theta,\rho,\sigma} \\ \mu u | \mathbf{y}=y \end{pmatrix} - Z_{X\Lambda} \begin{pmatrix} \hat{\beta}_{\theta,\rho,\sigma} \\ \mu u | \mathbf{y}=y \end{pmatrix} \right\|^2 \qquad (3.24)$$

$$= \left\| Y_{\Lambda} - Z_{X\Lambda} \begin{pmatrix} \hat{\beta}_{\theta,\rho,\sigma} \\ \mu u | \mathbf{y}=y \end{pmatrix} \right\|^2 + \left\| Z_{X\Lambda} \begin{pmatrix} \beta - \hat{\beta}_{\theta,\rho,\sigma} \\ u - \mu u | \mathbf{y}=y \end{pmatrix} \right\|^2 \qquad (3.25)$$

$$= \left\| Y_{\Lambda} - Z_{X\Lambda} \begin{pmatrix} \hat{\beta}_{\theta,\rho,\sigma} \\ \mu u | \mathbf{y}=y \end{pmatrix} \right\|^2 + \begin{pmatrix} \beta - \hat{\beta}_{\theta,\rho,\sigma} \\ u - \mu u | \mathbf{y}=y \end{pmatrix}^{\top} Z_{X\Lambda}^{\top} Z_{X\Lambda} \begin{pmatrix} \beta - \hat{\beta}_{\theta,\rho,\sigma} \\ u - \mu u | \mathbf{y}=y \end{pmatrix}. \qquad (3.26)$$

$Z_{X\Lambda}^{\top} Z_{X\Lambda}$ can be Cholesky decomposed as

$$Z_{X\Lambda}^{\top} Z_{X\Lambda} = \begin{pmatrix} R_X & \mathbf{0} \\ R_{ZX} & L_{\theta,\rho,\sigma}^{\top} \end{pmatrix}^{\top} \begin{pmatrix} R_X & \mathbf{0} \\ R_{ZX} & L_{\theta,\rho,\sigma}^{\top} \end{pmatrix}, \qquad (3.27)$$

where

$$L_{\theta,\rho,\sigma} L_{\theta,\rho,\sigma}^{\top} = Z_{\sigma\theta}^{\top} Z_{\sigma\theta} + \sqrt{\sigma_1^2 \sigma_2^2} \Sigma_u^{-1}. \qquad (3.28)$$

Thereafter,
$$
g(\beta, \mathbf{u}, \theta, \rho, \sigma) = \left\| Y_{\sigma} - X_{\sigma} \hat{\beta}_{\theta, \rho, \sigma} - Z_{\sigma\theta} \mu \mathbf{u} | \mathbf{y} = \mathbf{y} \right\|^2 + \sigma_1^2 \sigma_2^2 \mu_{\mathbf{u} | \mathbf{y} = \mathbf{y}}^{\top} \Sigma_{\mathbf{u}}^{-1} \mu_{\mathbf{u} | \mathbf{y} = \mathbf{y}} \\
+ \left\| R_X(\beta - \hat{\beta}_{\theta, \rho, \sigma}) \right\|^2 + \left\| R_{ZX}(\beta - \hat{\beta}_{\theta, \rho, \sigma}) + L_{\theta, \rho, \sigma}^{\top} (\mathbf{u} - \mu_{\mathbf{u} | \mathbf{y} = \mathbf{y}}) \right\|^2. \tag{3.29}
$$

By setting

$$
r(\hat{\beta}_{\theta, \rho, \sigma}, \mu_{\mathbf{u}|\mathbf{y}=\mathbf{y}}) = \left\| Y_{\sigma} - X_{\sigma} \hat{\beta}_{\theta, \rho, \sigma} - Z_{\sigma\theta} \mu_{\mathbf{u}|\mathbf{y}=\mathbf{y}} \right\|^2 + \sigma_1^2 \sigma_2^2 \mu_{\mathbf{u}|\mathbf{y}=\mathbf{y}}^{\top} \Sigma_{\mathbf{u}}^{-1} \mu_{\mathbf{u}|\mathbf{y}=\mathbf{y}}, \quad (3.30)
$$

and returning to the calculation of $f_{\mathbf{y}}(\mathbf{y})$, it comes

$$
\begin{align}
f_{\mathbf{y}}(\mathbf{y}) &= \frac{\int \exp \left[ - \frac{r(\hat{\beta}_{\theta,\rho,\sigma}, \mu_{\mathbf{u}|\mathbf{y}=\mathbf{y}}) + \|R_X(\beta-\hat{\beta}_{\theta,\rho,\sigma})\|^2 + \|R_{ZX}(\beta-\hat{\beta}_{\theta,\rho,\sigma}) + L_{\theta,\rho,\sigma}^{\top}(\mathbf{u}-\mu_{\mathbf{u}|\mathbf{y}=\mathbf{y}})\|^2}{2\sigma_1^2\sigma_2^2} \right] d\mathbf{u}}{(2\pi\sigma_1^2)^{N/2}(2\pi\sigma_2^2)^{N/2}(2\pi)^{q/2}|\Sigma_{\mathbf{u}}|^{1/2}} \nonumber \\
&= \frac{\exp \left[ - \frac{r(\hat{\beta}_{\theta,\rho,\sigma}, \mu_{\mathbf{u}|\mathbf{y}=\mathbf{y}}) + \|R_X(\beta-\hat{\beta}_{\theta,\rho,\sigma})\|^2}{2\sigma_1^2\sigma_2^2} \right]}{(2\pi\sigma_1^2)^{N/2}(2\pi\sigma_2^2)^{N/2}(2\pi)^{q/2}|\Sigma_{\mathbf{u}}|^{1/2}} \times \nonumber \\
&\qquad \int \exp \left[ - \frac{\|R_{ZX}(\beta-\hat{\beta}_{\theta,\rho,\sigma}) + L_{\theta,\rho,\sigma}^{\top}(\mathbf{u}-\mu_{\mathbf{u}|\mathbf{y}=\mathbf{y}})\|^2}{2\sigma_1^2\sigma_2^2} \right] d\mathbf{u}. \tag{3.31}
\end{align}
$$

By setting $v = R_{ZX}(\beta - \hat{\beta}_{\theta,\rho,\sigma}) + L_{\theta,\rho,\sigma}^{\top}(\mathbf{u} - \mu_{\mathbf{u}|\mathbf{y}=\mathbf{y}})$, $du = \frac{1}{|L_{\theta,\rho,\sigma}|}dv$ and

$$
\begin{align}
f_{\mathbf{y}}(\mathbf{y}) &= \frac{\exp\left[-\frac{r(\hat{\beta}_{\theta,\rho,\sigma}, \mu_{\mathbf{u}|\mathbf{y}=\mathbf{y}}) + \|R_X(\beta - \hat{\beta}_{\theta,\rho,\sigma})\|^2}{2\sigma_1^2\sigma_2^2}\right] (\sigma_1^2\sigma_2^2)^{\frac{q}{2}}}{(2\pi\sigma_1^2)^{N/2}(2\pi\sigma_2^2)^{N/2}|\Sigma_{\mathbf{u}}|^{1/2}|L_{\theta,\rho,\sigma}|} \nonumber \\
&\quad \times \int \frac{1}{(2\pi\sigma_1^2\sigma_2^2)^{\frac{q}{2}}} \exp\left[-\frac{\|\mathbf{v}\|^2}{2\sigma_1^2\sigma_2^2}\right] d\mathbf{v} \nonumber \\
&= \frac{\exp\left[-\frac{r(\hat{\beta}_{\theta,\rho,\sigma}, \mu_{\mathbf{u}|\mathbf{y}=\mathbf{y}}) + \|R_X(\beta - \hat{\beta}_{\theta,\rho,\sigma})\|^2}{2\sigma_1^2\sigma_2^2}\right] (\sigma_1^2\sigma_2^2)^{\frac{q}{2}}}{(2\pi\sigma_1^2)^{N/2}(2\pi\sigma_2^2)^{N/2}|\Sigma_{\mathbf{u}}|^{1/2}|L_{\theta,\rho,\sigma}|}. \tag{3.32}
\end{align}
$$

The log-likelihood to be maximized can therefore be expressed as,

$$
\begin{equation}
\begin{split}
\ell(\boldsymbol{\beta}, \boldsymbol{\theta}, \rho, \boldsymbol{\sigma} | \mathbf{y}) = & - \frac{r(\hat{\boldsymbol{\beta}}_{\boldsymbol{\theta},\rho,\sigma}, \mu_{\mathbf{u}|\mathbf{y}=\mathbf{y}}) + \|R_X(\boldsymbol{\beta} - \hat{\boldsymbol{\beta}}_{\boldsymbol{\theta},\rho,\sigma})\|^2}{2\sigma_1^2\sigma_2^2} \\
& - \frac{N-q}{2} \log(\sigma_1^2\sigma_2^2) \\
& - \frac{1}{2} \log(|\Sigma_{\mathbf{u}}|) - \frac{1}{2} \log(|L_{\boldsymbol{\theta},\rho,\sigma}|^2).
\end{split}
\tag{3.33}
\end{equation}
$$

□
By profiling out $\beta$, the partially profiled log-likelihood is

$$ \tilde{\ell}(\theta, \rho, \sigma | \mathbf{y}) = -\frac{r(\hat{\beta}_{\theta, \rho, \sigma}, \mu_{\mathbf{u}} | \mathbf{y}=\mathbf{y})}{2\sigma_1^2\sigma_2^2} - \frac{N-q}{2}\log(\sigma_1^2\sigma_2^2) \\ -\frac{1}{2}\log(|\Sigma_{\mathbf{u}}|) - \frac{1}{2}\log(|L_{\theta, \rho, \sigma}|^2), \quad (3.34) $$

replacing $\hat{\beta}_{\theta, \rho, \sigma}$ by $\beta$. Then, the partially profiled deviance comes

$$ -2\tilde{\ell}(\theta, \rho, \sigma | \mathbf{y}) = \frac{r(\hat{\beta}_{\theta, \rho, \sigma}, \mu_{\mathbf{u}} | \mathbf{y}=\mathbf{y})}{\sigma_1^2\sigma_2^2} + (N-q)\log(\sigma_1^2\sigma_2^2) \\ + \log(|\Sigma_{\mathbf{u}}|) + \log(|L_{\theta, \rho, \sigma}|^2). \quad (3.35) $$

This deviance is finally the criterion which will be minimized to obtaining the ML estimates of the parameters.

**Corollary 3.3.1.** Suppose that $\mathbf{y} = (y_1^\top, y_2^\top)^\top$ satisfies the bivariate linear mixed-effects model expressed by Equations (3.8 and 3.9). Taking into account the notations in the Theorem 4.2.1, the ML estimators $\hat{\beta}$, $\hat{\sigma}$, $\hat{\theta}$, $\hat{\rho}$ of $\beta$, $\sigma$, $\theta$ and $\rho$ satisfy

$$ (\hat{\theta}, \hat{\rho}, \hat{\sigma}) = \arg\max_{\theta, \rho, \sigma} \tilde{\ell}(\theta, \rho, \sigma | \mathbf{y}) \quad \text{and} \quad \hat{\beta} = \hat{\beta}_{\hat{\theta}, \hat{\rho}, \hat{\sigma}}. \quad (3.36) $$

### 3.3.2 REML criterion

By integrating the marginal density of $\mathcal{Y}$ with respect to the fixed effects, the REML criterion can be obtained [105]. This REML criterion is expressed through the following theorem.

**Theorem 3.3.2.** Suppose that $\mathbf{y} = (y_1^\top, y_2^\top)^\top$ satisfies the bivariate linear mixed-effects model expressed by Equations (3.8 and 3.9). Taking into account the notations in the Theorem 4.2.1, the REML criterion of $\sigma$, $\theta$ and $\rho$ given $\mathbf{y}$ is expressed as

$$ \mathcal{L}(\sigma, \theta, \rho | \mathbf{y}) = \frac{\exp\left[-\frac{r(\hat{\beta}_{\theta, \rho, \sigma}, \mu_{\mathbf{u}} | \mathbf{y}=\mathbf{y})}{2\sigma_1^2\sigma_2^2}\right] (\sigma_1^2\sigma_2^2)^{\frac{p+q-N}{2}}}{(2\pi)^{(2N-p)/2} |\Sigma_{\mathbf{u}}|^{1/2} |L_{\theta, \rho, \sigma}| |R_X|} \quad (3.37) $$

*Proof.*

$$ \mathcal{L}(\sigma, \theta, \rho | \mathbf{y}) = \int_{\mathbb{R}^p} f_{\mathcal{Y}}(\mathbf{y}) d\beta \quad (3.38) \\ = \frac{\exp\left[-\frac{r(\hat{\beta}_{\theta, \rho, \sigma}, \mu_{\mathbf{u}} | \mathbf{y}=\mathbf{y})}{2\sigma_1^2\sigma_2^2}\right] (\sigma_1^2\sigma_2^2)^{\frac{q}{2}}}{(2\pi\sigma_1^2)^{N/2} (2\pi\sigma_2^2)^{N/2} |\Sigma_{\mathbf{u}}|^{1/2} |L_{\theta, \rho, \sigma}|} \int_{\mathbb{R}^p} \exp\left[-\frac{\|R_X(\beta - \hat{\beta}_{\theta, \rho, \sigma})\|^2}{2\sigma_1^2\sigma_2^2}\right] d\beta $$
$$= \frac{\exp\left[-\frac{r(\hat{\beta}_{\theta, \rho, \sigma}, \mu_{\mathcal{U}|\mathbf{y}=\mathbf{y}})}{2\sigma_1^2\sigma_2^2}\right](\sigma_1^2\sigma_2^2)^{\frac{p+q-N}{2}}}{(2\pi)^{(2N-p)/2}|\Sigma_{\mathbf{u}}|^{1/2}|L_{\theta, \rho, \sigma}|R_X|} \int_{\mathbb{R}^p} (2\pi\sigma_1^2\sigma_2^2)^{-\frac{p}{2}} \exp\left[-\frac{\|\mathbf{t}\|^2}{2\sigma_1^2\sigma_2^2}\right] dt,$$

where $t = \beta - \hat{\beta}_{\theta, \rho, \sigma} \implies d\beta = \frac{1}{|R_X|}dt$ and $\int_{\mathbb{R}^p} (2\pi\sigma_1^2\sigma_2^2)^{-\frac{p}{2}} \exp\left[-\frac{\|\mathbf{t}\|^2}{2\sigma_1^2\sigma_2^2}\right] dt = 1$. $\square$

The REML criterion can also be expressed as

$$ \begin{aligned} \log(\mathcal{L}(\boldsymbol{\sigma}, \boldsymbol{\theta}, \rho|\mathbf{y})) &= \tilde{\mathcal{L}}(\boldsymbol{\sigma}, \boldsymbol{\theta}, \rho|\mathbf{y}) \\ &= -\frac{r(\hat{\beta}_{\theta, \rho, \sigma}, \mu_{\mathcal{U}|\mathbf{y}=\mathbf{y}})}{2\sigma_1^2\sigma_2^2} + \frac{p+q-N}{2}\log(\sigma_1^2\sigma_2^2) - \frac{1}{2}\log(|\Sigma_{\mathbf{u}}|) \\ &\quad -\frac{1}{2}\log(|L_{\theta, \rho, \sigma}|^2) - \frac{1}{2}\log(|R_X|^2), \end{aligned} \quad (3.39) $$

or as

$$ -2\tilde{\mathcal{L}}(\boldsymbol{\sigma}, \boldsymbol{\theta}, \rho|\mathbf{y}) = \frac{r(\hat{\beta}_{\theta, \rho, \sigma}, \mu_{\mathcal{U}|\mathbf{y}=\mathbf{y}})}{\sigma_1^2\sigma_2^2} + (N-p-q)\log(\sigma_1^2\sigma_2^2) + \log(|\Sigma_{\mathbf{u}}|) \\ + \log(|L_{\theta, \rho, \sigma}|^2) + \log(|R_X|^2), \quad (3.40) $$

which will be minimized to obtaining the REML estimates of the parameters.

## 3.4 Simulation studies

In this Section, the consistency of the estimators is proven through simulation studies, and we compare the present estimation procedure with the EM algorithm. For the sake of simplicity, these simulation studies are performed using simulated bivariate longitudinal data sets. In the following paragraph, we explain how we choose the parameters that have been used to simulate the 'working' data sets.

**The working data sets** We suppose that we are following up a sample of subjects where the goal is to evaluate how the growth of the weight and the height of the individuals of this population are jointly explained by the sex, the score of nutrition (Nscore) and the age. We randomly choose through a uniform distribution the score of nutrition between 20 and 50, and the age between 18 and 37, using the R software. All the computations in this paper are done using the R software. The subject's sex is also randomly chosen. The model under which we simulate the data sets is the following:

$n$ indicating the total number of subjects, for $i = 1, \dots, n$

$$ \begin{aligned} \text{weight}_i &= (1_{n_i}, \text{sex}_i, \text{Nscore}_i, \text{age}_i)\beta_1 + (1_{n_i}, \text{Nscore}_i)\gamma_{1i} + \varepsilon_{1i} \\ \text{height}_i &= (1_{n_i}, \text{sex}_i, \text{Nscore}_i, \text{age}_i)\beta_2 + (1_{n_i}, \text{Nscore}_i)\gamma_{2i} + \varepsilon_{2i} \end{aligned} \quad (3.41) $$
with

$$ \gamma_i = \begin{pmatrix} \gamma_{1i} \\ \gamma_{2i} \end{pmatrix} \sim \mathcal{N}(0, \bar{\Gamma}), \quad \varepsilon_{1i} \sim \mathcal{N}(0, \sigma_1^2 \mathbf{I}_{n_i}), \quad \varepsilon_{2i} \sim \mathcal{N}(0, \sigma_2^2 \mathbf{I}_{n_i}), \quad \gamma_i \perp \varepsilon_{1i} \perp \varepsilon_{2i} \quad (3.42) $$

The random effect related to the dependent variable 'weight' or 'height' is a vector composed by one random intercept and one random slope in the direction of the covariate 'Nscore'. The total number of observations is denoted by N.

We randomly choose $\beta_1, \beta_2, \sigma_1$ and $\sigma_2$ whose values are in the first column of Table 3.3. $\bar{\Gamma}$ is also randomly chosen such that it is positive definite, with the following form:

$$ \bar{\Gamma} = \begin{pmatrix} \eta_1^2 & \rho_{\eta\eta_1\eta_2} & \rho_{\eta_1\tau_1} & \rho_{\eta_1\tau_2} \\ \rho_{\eta\eta_1\eta_2} & \eta_2^2 & \rho_{\eta_2\tau_1} & \rho_{\eta_2\tau_2} \\ \rho_{\eta_1\tau_1} & \rho_{\eta_2\tau_1} & \tau_1^2 & \rho_{\tau\tau_1\tau_2} \\ \rho_{\eta_1\tau_2} & \rho_{\eta_2\tau_2} & \rho_{\tau\tau_1\tau_2} & \tau_2^2 \end{pmatrix} \quad (3.43) $$

In order to have an almost strong correlation between the marginal random effects, we set $\rho = 0.8$ and randomly choose all other parameters involved in the obtaining of $\bar{\Gamma}$. Thus, the obtained $\bar{\Gamma}$ is

$$ \bar{\Gamma} = \begin{pmatrix} 27.77 & 18.80 & 41.70 & 4.93 \\ 18.80 & 36.00 & 47.47 & 5.62 \\ 41.70 & 47.47 & 97.81 & 8.91 \\ 4.93 & 5.62 & 8.91 & 1.37 \end{pmatrix} \quad (3.44) $$

### 3.4.1 Estimates' performances

One practical way to show the consistency of an estimator is by computing its Mean Square Error (MSE). If the MSE of an estimator is asymptotically null, this estimator converges in probability, and is then consistent. In this Section, we gradually simulate data sets with larger sizes, $(N, n) \in \{(600, 50), (600, 60), \dots, (1000, 100), (1000, 300), \dots, (15000, 1000)\}$. We simulate one hundred data sets of each size and calculate the estimators' MSE using these data sets. This yields one hundred MSE for each size of dataset. This allows to compute the 95% CI (confidence interval) along with the mean of the MSE (one hundred mse) obtained for the hundred data sets of the same size. The results are contained in Table 3.1 and Table 3.2. The Table 3.1 shows that the asymptotique in the longitudinal data requires not only $n \to \infty$ and $N \to \infty$, but also $N/n \to \infty$. This means that it requires a sufficient total number of observations, a sufficient number of levels for the grouping factor and a sufficient number of observations for each level of the grouping factor. For example, in the Table 3.1, when the total number of observations is $N = 1000$, the MSE of $\bar{\Gamma}$ is better for $n = 100, 0.47 (0.03 - 1.29)$, than for $n = 300, 0.69 (0.06 - 1.94)$. Observing both Table 3.1 and Table 3.2 it is clear that, as the number of observations increase, the MSE descends to 0. We can conclude that the estimators constructed in this paper are consistent. The estimation procedure discussing in this paper may therefore be named Consistent estimates for the Multivariate Linear Mixed-Effects model (Cmlme). The Cmlme acronym will be used in the remainder of the paper for a question of simplicity.
<table><caption>Table 3.1: Mean Square Error of estimators with 95% CI estimated on 100 replications for values of $n \in \{50, 60, 100, 300\}$ and $N \in \{600, 1000, 3000\}$.</caption><thead><tr><th>Parameter</th><th>n</th><th>N = 600</th><th>N = 1000</th><th>N = 3000</th></tr></thead><tbody><tr><td rowspan="4">β<sub>1</sub></td><td>50</td><td>2.43 (0.11 - 7.11)</td><td>1.89 (0.22 - 4.89)</td><td>1.02 (0.07 - 2.44)</td></tr><tr><td>60</td><td>2.57 (0.14 - 7.87)</td><td>2.13 (0.26 - 5.55)</td><td>0.77 (0.10 - 2.09)</td></tr><tr><td>100</td><td>2.27 (0.16 - 5.61)</td><td>1.55 (0.14 - 5.17)</td><td>0.71 (0.04 - 1.85)</td></tr><tr><td>300</td><td>3.16 (0.14 - 10.54)</td><td>1.70 (0.10 - 4.55)</td><td>0.51 (0.04 - 1.36)</td></tr><tr><td rowspan="4">β<sub>2</sub></td><td>50</td><td>5.50(0.09 - 16.35)</td><td>3.26 (0.02 - 11.64)</td><td>2.06 (0.09 - 5.98)</td></tr><tr><td>60</td><td>5.06 (0.12 - 15.09)</td><td>3.22 (0.02 - 10.24)</td><td>1.78 (0.10 - 5.62)</td></tr><tr><td>100</td><td>4.33 (0.03 - 13.17)</td><td>2.37 (0.02 - 6.89)</td><td>1.06 (0.02 - 3.60)</td></tr><tr><td>300</td><td>4.58 (0.18 - 15.92)</td><td>2.43 (0.05 - 7.39)</td><td>0.90 (0.04 - 2.88)</td></tr><tr><td rowspan="4">σ<sub>1</sub></td><td>50</td><td>0.03 (0.00 - 0.14)</td><td>0.02 (0.00 - 0.09)</td><td>0.00 (0.00 - 0.02)</td></tr><tr><td>60</td><td>0.04 (0.00 - 0.11)</td><td>0.02 (0.00 - 0.08)</td><td>0.00 (0.00 - 0.02)</td></tr><tr><td>100</td><td>0.03 (0.00 - 0.13)</td><td>0.01 (0.00 - 0.06)</td><td>0.00 (0.00 - 0.01)</td></tr><tr><td>300</td><td>0.05 (0.00 - 0.23)</td><td>0.03 (0.00 - 0.13)</td><td>0.00 (0.00 - 0.02)</td></tr><tr><td rowspan="4">σ<sub>2</sub></td><td>50</td><td>0.04 (0.00 - 0.15)</td><td>0.03 (0.00 - 0.08)</td><td>0.00 (0.00 - 0.02)</td></tr><tr><td>60</td><td>0.04 (0.00 - 0.18)</td><td>0.03 (0.00 - 0.12)</td><td>0.01 (0.00 - 0.03)</td></tr><tr><td>100</td><td>0.06 (0.00 - 0.18)</td><td>0.03 (0.00 - 0.11)</td><td>0.01 (0.00 - 0.03)</td></tr><tr><td>300</td><td>0.13 (0.00 - 0.45)</td><td>0.04 (0.00 - 0.15)</td><td>0.01 (0.00 - 0.05)</td></tr><tr><td rowspan="4">Γ̄</td><td>50</td><td>0.90 (0.12 - 2.41)</td><td>0.62 (0.06 - 1.41)</td><td>0.45 (0.04 - 1.14)</td></tr><tr><td>60</td><td>1.07 (0.06 - 2.64)</td><td>0.68 (0.08 - 1.98)</td><td>0.25 (0.03 - 0.66)</td></tr><tr><td>100</td><td>0.90 (0.05 - 2.40)</td><td>0.47 (0.03 - 1.29)</td><td>0.21 (0.02 - 0.71)</td></tr><tr><td>300</td><td>1.45 (0.19 - 4.39)</td><td>0.69 (0.06 - 1.94)</td><td>0.17 (0.02 - 0.57)</td></tr></tbody></table>

<table><caption>Table 3.2: Mean Square Error of estimators with 95% CI estimated on 100 replications for values of $(n, N) \in \{(500, 7000), (600, 8000), (800, 10000), (1000, 15000)\}$.</caption><thead><tr><th>Parameter</th><th>n = 500, N = 7000</th><th>n = 600, N = 8000</th><th>n = 800, N = 10000</th><th>n = 1000, N = 15000</th></tr></thead><tbody><tr><td>β<sub>1</sub></td><td>0.22 (0.01 - 0.62)</td><td>0.17 (0.01 - 0.48)</td><td>0.16 (0.01 - 0.47)</td><td>0.11 (0.00 - 0.32)</td></tr><tr><td>β<sub>2</sub></td><td>0.32 (0.00 - 1.01)</td><td>0.34 (0.01 - 1.16)</td><td>0.25 (0.02 - 0.67)</td><td>0.21 (0.00 - 0.69)</td></tr><tr><td>σ<sub>1</sub></td><td>0.00 (0.00 - 0.00)</td><td>0.00 (0.00 - 0.01)</td><td>0.00 (0.00 - 0.00)</td><td>0.00 (0.00 - 0.00)</td></tr><tr><td>σ<sub>2</sub></td><td>0.00 (0.00 - 0.01)</td><td>0.00 (0.00 - 0.01)</td><td>0.00 (0.00 - 0.01)</td><td>0.00 (0.00 - 0.00)</td></tr><tr><td>Γ̄</td><td>0.09 (0.00 - 0.25)</td><td>0.07 (0.00 - 0.19)</td><td>0.06 (0.00 - 0.19)</td><td>0.03 (0.00 - 0.09)</td></tr></tbody></table>
3.4.2 Comparison with EM-based estimates

In this Section, we compare the estimation procedure based on EM algorithm with the Cmlme. This
comparison is performed regarding the accuracy of the estimates, whether or not the starting values
of the two algorithms (EM and Cmlme) are naive or advised. We mean by naive starting values,
those which are randomly chosen (without specific control), and by advised starting values those
obtained by fitting separately each dimension of the bivariate model. The results of these marginal
fitting are, indeed, the advised starting values for the bivariate model estimation procedure. The
starting values are the same for both Cmlme and EM algorithms. The number of iteration required
for convergence, for each algorithm, is also discussed. Our methodology consists in simulating
thirty longitudinal data sets of size (N = 3000, n = 300) and fit the model to each of these data
sets using the EM algorithm and the Cmlme, respectively. This allows to compute both the 95%
CI and the empirical mean of the thirty estimates in each case (naive and advised starting values).
The obtained results are in Table 3.3 and Table 3.4. Table 3.3 contains the empirical means of
the estimates with their 95% CI, and the minimum, the maximum and the average number of
iterations. Table 3.4 contains the empirical relative error of the estimators with their 95% CI.
These results show that in the case of naive initialization, the Cmlme estimators outperform the
EM estimators. For example, the component of β₁ which is 14.00 is well estimated by Cmlme,
14.02 (13.27 – 14.45) with an empirical relative error of 0.02 (0.00 – 0.04), but poorly estimated
by EM, -2.05 (-4.70 – -0.40) with an empirical relative error of 1.14 (1.01 – 1.32). In the
case of advised initialization, both Cmlme and EM algorithms perform well, but Cmlme converge
faster (64 iterations in average) than EM (169 iterations in average). The number of iteration
required by Cmlme with advised initializations range from 48 to 89 and from 56 to 103 for naive
initializations. The Cmlme with naive initialization therefore needs more iterations than Cmlme
with advised initialization for converge. This is expected and may be explained by the fact that
the advised initialization values contain some information from the data of interest, and the naive
starting points do not.

The empirical mean of the random effects covariance matrix, $\bar{\Gamma}$, is well estimated with advised initializations:

$$
\bar{\Gamma}_{\text{adv}}^{\text{Cmlme}} = \begin{pmatrix} 25.74 & 16.34 & 35.75 & 4.50 \\ 16.34 & 34.83 & 43.97 & 5.37 \\ 35.75 & 43.97 & 82.44 & 8.43 \\ 4.50 & 5.37 & 8.43 & 1.32 \end{pmatrix}, \quad \text{with} \quad \sigma_{\Gamma_{\text{adv}}}^{\text{Cmlme}} = \begin{pmatrix} 5.73 & 2.74 & 4.45 & 0.61 \\ 2.74 & 2.59 & 3.96 & 0.44 \\ 4.45 & 3.96 & 12.07 & 0.80 \\ 0.61 & 0.44 & 0.80 & 0.11 \end{pmatrix} \tag{3.45}
$$

and

$$
\bar{\Gamma}_{\text{adv}}^{\text{EM}} = \begin{pmatrix} 23.66 & 16.45 & 32.60 & 4.45 \\ 16.45 & 34.82 & 43.39 & 5.39 \\ 32.61 & 43.39 & 75.16 & 8.66 \\ 4.45 & 5.39 & 8.65 & 1.31 \end{pmatrix}, \quad \text{with} \quad \sigma_{\Gamma_{\text{adv}}}^{\text{EM}} = \begin{pmatrix} 7.49 & 2.71 & 5.29 & 0.76 \\ 2.71 & 2.59 & 3.71 & 0.45 \\ 5.30 & 3.71 & 13.36 & 0.80 \\ 0.76 & 0.45 & 0.80 & 0.11 \end{pmatrix} \tag{3.46}
$$
<table><caption>Table 3.3: EM estimates compared with Cmlme estimates on the same data sets. Empirical estimates with their 95% CI and the number of iteration required for convergence.</caption><thead><tr><th rowspan="3">Parameter</th><th rowspan="3">Value</th><th colspan="6">Naive initialization</th><th colspan="2">Advised initialization</th></tr><tr><th colspan="2">Cmlme</th><th colspan="2">EM</th><th colspan="2">Cmlme</th><th colspan="2">EM</th></tr><tr><th>Emp. Mean</th><th>95% CI</th><th>Emp. Mean</th><th>95% CI</th><th>Emp. Mean</th><th>95% CI</th><th>Emp. Mean</th><th>95% CI</th></tr></thead><tbody><tr><td rowspan="3">β1</td><td>50.67</td><td>50.79</td><td>49.14 - 52.11</td><td>13.47</td><td>-76.70 - 43.37</td><td>50.80</td><td>49.15 - 52.12</td><td>50.78</td><td>49.09 - 52.01</td></tr><tr><td>-4.80</td><td>-5.00</td><td>-8.39 - -3.66</td><td>-4.79</td><td>-8.08 - -3.42</td><td>-5.02</td><td>-8.39 - -3.66</td><td>-4.98</td><td>-8.38 - -3.65</td></tr><tr><td>14.00</td><td>14.02</td><td>13.27 - 14.45</td><td>2.05</td><td>-4.70 - -0.40</td><td>14.02</td><td>13.28 - 14.45</td><td>14.02</td><td>13.27 - 14.45</td></tr><tr><td rowspan="3">β2</td><td>2.70</td><td>2.70</td><td>2.66 - 2.72</td><td>2.69</td><td>2.66 - 2.72</td><td>2.70</td><td>2.66 - 2.72</td><td>2.70</td><td>2.66 - 2.72</td></tr><tr><td>13.20</td><td>13.65</td><td>11.79 - 15.06</td><td>-84.47</td><td>-114.28 - -50.63</td><td>13.65</td><td>11.79 - 15.07</td><td>13.68</td><td>11.81 - 15.14</td></tr><tr><td>-2.80</td><td>-2.80</td><td>-4.74 - -0.43</td><td>-2.75</td><td>-4.90 - 0.21</td><td>-2.81</td><td>-4.79 - -0.43</td><td>-2.85</td><td>-4.80 - -0.51</td></tr><tr><td rowspan="3">σ1</td><td>27.00</td><td>27.00</td><td>26.87 - 27.10</td><td>0.90</td><td>-1.62 - 2.68</td><td>27.00</td><td>26.87 - 27.10</td><td>27.00</td><td>26.87 - 27.10</td></tr><tr><td>1.70</td><td>1.68</td><td>1.64 - 1.71</td><td>1.68</td><td>1.64 - 1.71</td><td>1.68</td><td>1.64 - 1.71</td><td>1.68</td><td>1.64 - 1.71</td></tr><tr><td>5.80</td><td>5.79</td><td>5.62 - 5.92</td><td>5.78</td><td>5.64 - 5.98</td><td>5.78</td><td>5.64 - 5.92</td><td>5.79</td><td>5.65 - 5.94</td></tr><tr><td rowspan="3">σ2</td><td>7.60</td><td>7.61</td><td>7.34 - 7.74</td><td>7.59</td><td>7.33 - 7.73</td><td>7.61</td><td>7.34 - 7.73</td><td>7.63</td><td>7.36 - 7.73</td></tr><tr><td colspan="2">Min</td><td>56</td><td>63</td><td colspan="2">-</td><td>48</td><td>-</td><td>14</td><td>-</td></tr><tr><td colspan="2">Mean</td><td>71</td><td>109</td><td colspan="2">-</td><td>64</td><td>-</td><td>169</td><td>-</td></tr><tr><td colspan="2">Max</td><td>103</td><td>157</td><td colspan="2">-</td><td>89</td><td>-</td><td>645</td><td>-</td></tr></tbody></table>

$\sigma_{\Gamma_{\text{adv}}}^{\text{Cmlme}}$ and $\sigma_{\Gamma_{\text{adv}}}^{\text{EM}}$ contain the standard deviations of the entries of $\bar{\Gamma}_{\text{adv}}^{\text{Cmlme}}$ and $\bar{\Gamma}_{\text{adv}}^{\text{EM}}$, respectively. It seems that the empirical standard deviation of the higher entries of $\bar{\Gamma}$ are bigger with EM than with Cmlme. For example, the standard deviation of $\bar{\Gamma}_{11} = 27.77$ is 7.49 for EM, but 5.73 for Cmlme. Same remark about the standard deviations of $\bar{\Gamma}_{31}$, $\bar{\Gamma}_{32}$ and $\bar{\Gamma}_{33}$, comparing EM and Cmlme. This may be explained by the fact that Cmlme estimators are more consistent than EM
<table><caption>Table 3.4: EM estimates compared with Cmlme estimates on the same data sets. Empirical relative error of the estimates with their 95% CI</caption><thead><tr><th rowspan="3">Parameter</th><th rowspan="3">Value</th><th colspan="4">Naive initialization</th><th colspan="4">Advised initialization</th></tr><tr><th colspan="2">Cmlme</th><th colspan="2">EM</th><th colspan="2">Cmlme</th><th colspan="2">EM</th></tr><tr><th>R. Error</th><th>95% CI</th><th>R. Error</th><th>95% CI</th><th>R. Error</th><th>95% CI</th><th>R. Error</th><th>95% CI</th></tr></thead><tbody><tr><td rowspan="4">β<sub>1</sub></td><td>50.67</td><td>0.01</td><td>0.00 - 0.03</td><td>0.73</td><td>0.01 - 2.18</td><td>0.01</td><td>0.00 - 0.03</td><td>0.01</td><td>0.00 - 0.03</td></tr><tr><td>-4.80</td><td>0.21</td><td>0.02 - 0.32</td><td>0.21</td><td>0.00 - 0.32</td><td>0.21</td><td>0.02 - 0.32</td><td>0.21</td><td>0.00 - 0.33</td></tr><tr><td>14.00</td><td>0.02</td><td>0.00 - 0.04</td><td>1.14</td><td>1.01 - 1.32</td><td>0.02</td><td>0.00 - 0.04</td><td>0.02</td><td>0.00 - 0.04</td></tr><tr><td>2.70</td><td>0.00</td><td>0.00 - 0.01</td><td>0.00</td><td>0.00 - 0.01</td><td>0.00</td><td>0.00 - 0.01</td><td>0.00</td><td>0.00 - 0.01</td></tr><tr><td rowspan="4">β<sub>2</sub></td><td>13.20</td><td>0.07</td><td>0.00 - 0.14</td><td>7.39</td><td>4.19 - 9.40</td><td>0.07</td><td>0.00 - 0.14</td><td>0.07</td><td>0.00 - 0.14</td></tr><tr><td>-2.80</td><td>0.43</td><td>0.00 - 0.84</td><td>0.43</td><td>0.02 - 1.07</td><td>0.43</td><td>0.00 - 0.84</td><td>0.43</td><td>0.00 - 0.81</td></tr><tr><td>27.00</td><td>0.00</td><td>0.00 - 0.00</td><td>0.96</td><td>0.88 - 1.02</td><td>0.00</td><td>0.00 - 0.00</td><td>0.00</td><td>0.00 - 0.00</td></tr><tr><td>1.70</td><td>0.01</td><td>0.00 - 0.02</td><td>0.01</td><td>0.00 - 0.02</td><td>0.01</td><td>0.00 - 0.02</td><td>0.01</td><td>0.00 - 0.03</td></tr><tr><td>σ<sub>1</sub></td><td>5.80</td><td>0.01</td><td>0.00 - 0.02</td><td>0.01</td><td>0.00 - 0.03</td><td>0.01</td><td>0.00 - 0.02</td><td>0.01</td><td>0.00 - 0.03</td></tr><tr><td>σ<sub>2</sub></td><td>7.60</td><td>0.01</td><td>0.00 - 0.02</td><td>0.01</td><td>0.00 - 0.04</td><td>0.01</td><td>0.00 - 0.02</td><td>0.01</td><td>0.00 - 0.02</td></tr></tbody></table>

estimators. In the case of naive initializations, Cmlme provides a well estimated empirical mean of $\bar{\Gamma}$, when the estimated $\bar{\Gamma}$ provided by EM is very bad (we choose not to show it here).
$$ \bar{\Gamma}_{\text{naiv}}^{\text{Cmlme}} = \begin{pmatrix} 24.76 & 16.18 & 34.43 & 4.47 \\ 16.18 & 34.77 & 43.91 & 5.35 \\ 34.43 & 43.91 & 80.67 & 8.40 \\ 4.47 & 5.35 & 8.40 & 1.32 \end{pmatrix}, \quad \text{with} \quad \sigma_{\bar{\Gamma}_{\text{naiv}}^{\text{Cmlme}}} = \begin{pmatrix} 7.58 & 2.67 & 6.76 & 0.61 \\ 2.67 & 2.54 & 3.91 & 0.42 \\ 6.76 & 3.91 & 13.37 & 0.80 \\ 0.61 & 0.42 & 0.80 & 0.11 \end{pmatrix} \quad (3.47) $$

$\bar{\Gamma}_{\text{naiv}}^{\text{Cmlme}}$ compared to $\bar{\Gamma}_{\text{adv}}^{\text{Cmlme}}$ and $\sigma_{\bar{\Gamma}_{\text{naiv}}^{\text{Cmlme}}}$ compared to $\sigma_{\bar{\Gamma}_{\text{adv}}^{\text{Cmlme}}}$ show a slight difference which reveals a tiny sensibility of Cmlme to the starting values. This may be corrected by doing more than one evaluation of the model's deviance.

For all the simulation studies, we use the ML deviance criterion (Equation 3.35) and have minimized it using the nlminb function under R software. Thus, the estimates obtained are from the ML estimators. In this paper, we do not provide an application of REML estimates.

## 3.5 Application on malaria dataset

### 3.5.1 Data description

The data that we analyze here come from a study which was conducted in 9 villages (Avamé centre, Gbédjougo, Houngo, Anavié, Dohinoko, Gbétaga, Tori Cada Centre, Zébè and Zoungoudo) of Tori Bossito area (Southern Benin), where *P. falciparum* is the commonest species in the study area (95%) [47] from June 2007 to January 2010. The aim of this study was to evaluate the determinants of malaria incidence in the first months of life of child in Benin.

Mothers ($n = 620$) were enrolled at delivery and their newborns were actively followed-up during the first year of life. One questionnaire was conducted to gather information on women's characteristics (age, parity, use of Intermittent Preventive Treatment during pregnancy (IPTp) and bed net possession) and on the course of their current pregnancy. After delivery, thick and thin placental blood smears were examined to detect placental infection defined by the presence of asexual forms of *P. falciparum*. Maternal peripheral blood as well as cord blood were collected. At birth, newborn's weight and length were measured and gestational age was estimated.

During the follow-up of newborns, axillary temperature was measured weekly. In case of temperature higher than 37.5°C, mothers were told to bring their children to the health center where a questionnaire was filled out. A rapid diagnostic test (RDT) for malaria was performed and a thick blood smear (TBS) made. Symptomatic malaria cases, defined as fever (> 37.5°C) with TBS and/or RDT positive, were treated with an artemisinin-based combination. Systematically, TBS were made every month to detect asymptomatic infections. Every three months, venous blood was sampled to quantify the level of antibody against malaria promised candidate vaccine antigens. Finally, the environmental risk of exposure to malaria was modeled for each child, derived from a statistical predictive model based on climatic, entomological parameters, and characteristics of children's immediate surroundings. Also every 3 months (at 3, 6, 9, 12, 15, 18 months 130 of age), infant blood samples were collected.

Concerning the antibody quantification, two recombinant *P. falciparum* antigens where used to perform IgG subclass (IgG1 and IgG3) antibody. Recombinants antigens MSP2 (3D7 and
FC27) were from La Trobe University [7, 124]. GLURP-R0 (amino acids 25-514, F32 strain) and GLURP-R2 (amino acids 706-1178, 140 F32 strain) were also expressed. The antibodies were quantified in plasma at different times and ADAMSEL FLPb039 software (http://www.malariaresearch.eu/content/software) was used to analyze automatically the ELISA optical density (OD) leading to antibody concentrations in (µg/mL).

In this paper, we use some of the data and we rename the proteins used in the study, for reasons of the protection of these data. Thus, the proteins we use here, are named A1, A2, B and C, and are related to the antigens IgG1 and IgG3 as mentioned above. Information contained in the multivariate longitudinal dataset of malaria are described in the Table 3.5, where Y denotes an antigen which is one of the following:

$$ \text{IgG1\_A1, IgG3\_A1, IgG1\_A2, IgG3\_A2, IgG1\_B, IgG3\_B, IgG1\_C, IgG3\_C} \quad (3.48) $$

<table><caption>Table 3.5: Variables present in the empirical dataset</caption><thead><tr><th>N°</th><th>Variable</th><th>Description</th></tr></thead><tbody><tr><td>1</td><td>id</td><td>Child ID</td></tr><tr><td>2</td><td>conc.Y</td><td>concentration of Y</td></tr><tr><td>3</td><td>conc_CO.Y</td><td>Measured concentration of Y in the umbilical cord blood</td></tr><tr><td>4</td><td>conc_M3.Y</td><td>Predicted concentration of Y in the child's peripheral blood at 3 months</td></tr><tr><td>5</td><td>ap</td><td>Placental apposition</td></tr><tr><td>6</td><td>hb</td><td>Hemoglobin level</td></tr><tr><td>7</td><td>inf_trim</td><td>Number of malaria infections in the previous 3 months</td></tr><tr><td>8</td><td>pred_trim</td><td>Quarterly average number of mosquitoes child is exposed to</td></tr><tr><td>9</td><td>nutri_trim</td><td>Quarterly average nutrition scores</td></tr></tbody></table>

### 3.5.2 Data analysis

The aim of the analysis of these data is to evaluate the effect of the malaria infection on the child's immune acquisition (against malaria). Since the antigens which characterize the child's immune status interact together in the human body, we analyze the characteristics of the joint distribution of these antigens, conditionally to the malaria infection and other factors of interest. The dependent variables are then provided by conc.Y (Table 3.5) which describes the level of the antigen Y in the children at 3, 6, 9, 12, 15 and 18 months. All other variables in the Table 3.5 are covariates. We then have 8 dependent variables which describe the longitudinal profile (in the child) of the proteins listed in Equation 3.48.

To illustrate the stability of our approach, we are fitting here a bivariate model to the data, with IgG1\_A1 and IgG3\_A2 as dependent variables:

$$ \text{conc.IgG1\_A1} = (1, ap, \text{conc\_CO.IgG1\_A1}, \text{conc\_M3.IgG1\_A1}, \text{hb}, \text{inf\_trim}, \text{pred\_trim}, \text{nutri\_trim})\beta_1 + (1, \text{inf\_trim})\gamma_1 + \varepsilon_1 $$
$$ \text{conc.IgG3\_A2} = (\mathbf{1}, \text{ap}, \text{conc\_CO.IgG3\_A2}, \text{conc\_M3.IgG3\_A2}, \text{hb}, \text{inf\_trim}, \\ \text{pred\_trim}, \text{nutri\_trim})\beta_2 + (\mathbf{1}, \text{inf\_trim})\gamma_2 + \varepsilon_2, \quad (3.49) $$

with

$$ \gamma = (\gamma_1^\top, \gamma_2^\top)^\top \sim \mathcal{N}(\mathbf{0}, \bar{\Gamma}), \quad \varepsilon = (\varepsilon_1^\top, \varepsilon_2^\top)^\top \sim \mathcal{N}(\mathbf{0}, \begin{pmatrix} \sigma_1^2\mathbf{I} & 0 \\ 0 & \sigma_2^2\mathbf{I} \end{pmatrix}). \quad (3.50) $$

Our strategy is to 1) fit the model to the data by running the Cmlme algorithm using 25 different naive starting points and 2) retain the estimates related to the best likelihood (the minimum of the 25 deviances) as the true parameters and compute the estimators' MSE using the 24 others estimates. This may allow to evaluate how much the Cmlme algorithm is sensitive to the starting points. The results are contained in the Table 3.6.

<table><caption>Table 3.6: Empirical data analysis.</caption><thead><tr><th rowspan="3">Covariates</th><th colspan="4">Response variables</th></tr><tr><th colspan="2">conc.IgG1_A1</th><th colspan="2">conc.IgG3_A2</th></tr><tr><th>Estimate</th><th>MSE</th><th>Estimate</th><th>MSE</th></tr></thead><tbody><tr><td>Intercept</td><td>0.609</td><td>9.05 × 10<sup>-5</sup></td><td>-1.626</td><td>3.18 × 10<sup>-5</sup></td></tr><tr><td>ap</td><td>-0.093</td><td>1.06 × 10<sup>-5</sup></td><td>-0.337</td><td>1.04 × 10<sup>-6</sup></td></tr><tr><td>conc_CO.IgG1_A1</td><td>0.160</td><td>1.68 × 10<sup>-6</sup></td><td>—</td><td>—</td></tr><tr><td>conc_M3.IgG1_A1</td><td>0.148</td><td>9.85 × 10<sup>-6</sup></td><td>—</td><td>—</td></tr><tr><td>conc_CO.IgG3_A2</td><td>—</td><td>—</td><td>0.047</td><td>6.44 × 10<sup>-7</sup></td></tr><tr><td>conc_M3.IgG3_A2</td><td>—</td><td>—</td><td>0.155</td><td>2.22 × 10<sup>-7</sup></td></tr><tr><td>hb</td><td>-0.162</td><td>3.22 × 10<sup>-7</sup></td><td>-0.345</td><td>1.35 × 10<sup>-7</sup></td></tr><tr><td>inf\_trim</td><td>0.369</td><td>1.89 × 10<sup>-6</sup></td><td>0.696</td><td>5.09 × 10<sup>-7</sup></td></tr><tr><td>pred\_trim</td><td>-0.003</td><td>5.25 × 10<sup>-8</sup></td><td>0.017</td><td>1.49 × 10<sup>-8</sup></td></tr><tr><td>nutri\_trim</td><td>0.024</td><td>5.81 × 10<sup>-6</sup></td><td>0.115</td><td>3.71 × 10<sup>-5</sup></td></tr><tr><td>σ<sub>1</sub> and σ<sub>2</sub></td><td>1.395</td><td>4.96 × 10<sup>-6</sup></td><td>1.626</td><td>2.42 × 10<sup>-5</sup></td></tr></tbody></table>

Based on these results, the influence of the starting points on the Cmlme algorithm is very low (see the MSE in Table 3.6). The estimated random effects covariance matrix is

$$ \Gamma = \begin{pmatrix} 0.58 & -0.13 & 0.74 & -0.36 \\ -0.13 & 0.23 & -0.39 & 0.37 \\ 0.74 & -0.39 & 0.94 & -0.24 \\ -0.36 & 0.37 & -0.24 & 0.34 \end{pmatrix} \quad (3.51) $$

with an MSE of 0.0095.

## 3.6 Conclusion

In the context of fitting multivariate linear mixed-effects model having homoscedastic dimensional residuals, we have suggested ML and REML estimation strategies by profiling the model's
deviance and Cholesky factorizing the random effect covariance matrix. This approach can be considered as the generalization of the approach used by [16] in the R software lme4 package. Through extensive simulation studies, we have illustrated that the present approach outperforms the traditional EM estimates and provides estimates that are consistent for both fixed effects and variance components. Another interesting characteristic is its robustness relative to the initial value of the optimization procedure which can be randomly chosen without affecting the estimation results. Furthermore, the profiled ML or REML criterion's optimization can be easily and rapidly performed using an existing optimizer in the R software. Further considerations of this approach may include heteroscedastic residuals as well as residuals correlated with random effects, where the theoretical consistency of the resulting estimators will be demonstrated.
# FIXED EFFECTS SELECTION IN THE LINEAR MIXED-EFFECTS MODEL USING ADAPTIVE RIDGE PROCEDURE FOR L0 PENALTY PERFORMANCE

Dans ce chapitre, nous faisons de la sélection d'effets fixes dans le modèle linéaire à effets mixtes (unidimensionnel). L'approche utilisée est une optimisation de la déviance pénalisée, où le terme de pénalité est de type $L_2$ en présence d'une matrice de poids itérativement mise à jour. C'est ce que nous avons appelé la procédure *Adaptive Ridge itérative* dont la qualité de selection approxime celle résultant d'une pénalité de type $L_0$. Ici, le calcul de la vraisemblance a été inspiré de celui proposé par Bates et al. (2014) avec une modification qui permet de "libérer" le vecteur des effets fixes afin de permettre une sélection plus efficace de ces effets. Les résultats obtenus sur des données simulées sont de très bonnes qualités, aussi bien en terme de sélection que d'estimation. Le présent chapitre fait aussi l'objet d'un article de recherche soumis dans une revue scientifique.
Fixed effects selection in the linear mixed-effects model using adaptive
ridge procedure for L₀ penalty performance

Eric Hougla Adjakossa¹,²*, Mahouton Norbert Hounkonnou², Gregory Nuel¹

1 Laboratoire de Probabilités et Modèles Aléatoires /Université Pierre et Marie Curie, Case courrier 188 - 4, Place Jussieu 75252 Paris cedex 05 France

2 University of Abomey-Calavi, 072 B.P. 50 Cotonou, Republic of Benin

* ericadjakossah@gmail.com

abstract

This paper is concerned with the selection of fixed effects along with the estimation of fixed ef-
fects, random effects and variance components in the linear mixed-effects model. We introduce
a selection procedure based on an adaptive ridge (AR) penalty of the profiled likelihood, where
the covariance matrix of the random effects is Cholesky factorized. This selection procedure is in-
tended to both low and high-dimensional settings where the number of fixed effects is allowed to
grow exponentially with the total sample size, yielding technical difficulties due to the non-convex
optimization problem induced by L₀ penalties. Through extensive simulation studies, the proce-
dure is compared to the LASSO selection and appears to enjoy the model selection consistency as
well as the estimation consistency.

**4.1 Introduction**

During the last two decades, selection procedures in the linear mixed-effects model have been an
active research topic due to the appealing features of the model and the advent of modern technolo-
gies facilitating the collection of many variables in scientific studies. Many of these variables are
typically included in the full model at the initial stage of modeling to reduce model approximation
error, and due to the complexity of the mixed-effects models, inferences and interpretations of the
estimated models become challenging as the dimension of fixed or random effects increases [54].
The selection of important fixed or random effects has thus become a fundamental problem in the
analysis of grouped data using mixed-effects models, especially in the high-dimensional settings
where the fixed or the random effects vector dimension is allowed to grow exponentially with the
sample size.

Generally, model selection procedures can be viewed as covering three main approaches: the
hypothesis testing procedures, the regularization procedures and other procedures which include
the Bayesian selection methods. The testing procedures include the ordinary hypothesis tests and
the selection methods based on generalized information criteria. Examples of using testing hy-
pothesis for models selection in the mixed-effects model context include Lin’s works [112] who
proposed a simple global variance component tests, which are locally asymptotically most pre-
cise and are robust in the sense that no assumption about the parametric form of the random
effects is made. Despite their global form expressions which require only the fitting of conven-
tional generalized linear models, the critical values of the global test statistics, which are based
on large sample theory, are less accurate when the number of levels of each random effect is
small, e.g. less than 15. Edwards and his co-workers [50] extended the traditional coefficient of
determination R² for the linear mixed-effects model y ~ N(Xβ, Σ = ZΓZ' + σ²I_N), where
they introduced a statistic R<sub>β</sub><sup>2</sup> = (q - 1)ν<sup>-1</sup>F(β̂, Σ̂)/[1 + (q - 1)ν<sup>-1</sup>F(β̂, Σ̂)], with F(β̂, Σ̂) =
(Cβ)<sup>⊤</sup> [C(X'Σ<sup>-1</sup>X)<sup>-1</sup>C<sup>⊤</sup>]<sup>-1</sup> Cβ/rank(C), ν = N - rank(X) = N - q, C = [0<sub>(q-1)×1</sub>I<sub>q-1</sub>]
of rank q - 1, in testing H<sub>0</sub> : Cβ = 0. R<sub>β</sub><sup>2</sup> measures the multivariate association between the
repeated outcomes and the fixed effects in the context of longitudinal data analysis. This R<sub>β</sub><sup>2</sup> statis-
tic arises as a 1 - 1 function of an appropriate F statistic (i.e., F(β̂, Σ̂)) for testing all the fixed
effects, except the intercept. More precisely, R<sub>β</sub><sup>2</sup> compares the full model with a null model having
no fixed effect except typically the intercept. R<sub>β</sub><sup>2</sup> is then generalized to define a partial R<sup>2</sup> statistic
for marginal fixed effects of all sorts. Although this testing-based selection procedure of fixed ef-
fects is very useful, one of its major drawback is that the choice of the denominator of R<sub>β</sub><sup>2</sup> clearly
affects the rate of convergence as N → ∞, and may change the parameter being estimated. Exam-
ples of R<sup>2</sup>-based selection of fixed effects in the mixed-effects model include [178], [202], [101]
and references therein, where generally, too much restrictions are made on the random effects
covariance matrix, and clearly may not be appropriate for a wide range of data analysis.

Based on testing procedures, a stepwise procedure can be constructed for selecting important
fixed or random effects using generalized information criteria [54, 144], which are a generalization
of Akaike's information criterion (AIC) [4] and the Bayesian information criterion (BIC) [167].
[144] extended the Generalized information criterion (GIC) proposed by [150] in order to con-
struct a procedure for selecting fixed and random effects in the linear mixed-effects model. Here,
following [173], the asymptotic behavior of the extended GIC method for selecting fixed effects is
studied, and the results from simulations show that if the signal-to-noise ratio is moderate or high,
the percentages of choosing the correct fixed effects by the GIC procedure are close to one for fi-
nite samples. Another examples of GIC like selection procedures include those proposed by [122],
[85] and [21]. These strategies suggest a unified approach for choosing a parameters vector β that
maximizes the penalized likelihood

$$
n^{-1} \ell_n(\beta|y) - \sum_{j=1}^{p} p_{\lambda}(|\beta_j|) \quad (4.1)
$$

which arose from the Kullback-Leibler (KL) divergence $-\ell_n(\hat{\beta}|y) + \lambda\|\beta\|_0$ of the fitted model
from the true model [3], where $\ell_n(\cdot|y)$ is the log-likelihood function, $\hat{\beta}$ is the maximum likelihood
estimator of $\beta$ and $p_\lambda$ is the penalty function indexed by the regularization parameter $\lambda \ge 0$. The
$L_0$-norm $\|\beta\|_0$ of $\beta$ counts the number of non-vanishing components ($\beta_j \ne 0$) in $\beta$, and arises
naturally in many classical model selection methods. Although involving a nice interpretation of
the best subset selection, and admitting good sampling properties [11], its computation is infeasible in high dimensional statistical endeavors [53]. Other penalty functions are regularly used in the literature. A natural generalization of $L_0$ penalty is the so-called bridge (noted $L_q$) penalty in [63], where $p_\lambda(|\beta|) = \lambda|\beta||_{L_q}^q$ for $0 < q \le 2$. The $L_q$ penalty encompasses $L_0$, $L_1$ - LASSO [186] - and $L_2$ (ridge) penalties. Since none of the $L_q$ penalties satisfies all the required properties (sparsity, approximate unbiasedness and continuity, see [52] for more details) for their resulting parameter estimators, other penalties including SCAD [51, 52] and MCP [207] are introduced in the literature.

Testing-based stepwise selection procedures, where $\lambda$ is fixed ($\lambda = 1$ for AIC, $\lambda = \log(n)/2$ for BIC, $\lambda = \log(\log n)/2$ for HQIC [85], $\lambda = (\log n + 1)/2$ for CAIC [21], for example) are computationally expensive for high-dimensional settings and ignore stochastic errors in the variable selection process [52]. Another severe drawback is their lack of stability [24]. Penalized likelihood approaches using data-driven choice of $\lambda$ are generally preferred to handle the high-dimensional selection problem. These methods are referred to as regularization procedures. Here, the regularization parameter $\lambda$ (also called the tuning parameter) is generally estimated using cross-validation methods [23, 52, 66, 186]. [95] consider the selection of both fixed and random effects in a general class of mixed effects models by optimizing the penalized likelihood criterion $Q_\lambda(\theta|\theta_{old}) = Q(\theta|\theta_{old}) - n \sum_{j=1}^p p_\lambda(|\beta_j|)$ using the EM algorithm [44], where $\theta$ is the vector of all the unknown parameters, $Q(\cdot|\theta_{old})$ is the resulting function of the EM algorithm E-step, and $p_\lambda$ is either the SCAD or the adaptive lasso - ALASSO - [209] penalty. They approximate the integral in $Q(\theta|\theta_{old})$ by using a Markov chain Monte Carlo and introduce the $IC_Q(\lambda) = -2Q(\hat{\theta}_\lambda|\hat{\theta}_0) + c_n(\hat{\theta}_\lambda)$ statistic [94] for selecting the regularization parameter, where $\hat{\theta}_0 = \arg\max_\theta \ell(\theta)$ is the unpenalized maximum likelihood estimate and $c_n(\theta)$ a function of the data and fitted model.

[54] introduce a class of variable selection methods for fixed effects using a penalized profiled likelihood method, where the unknown covariance matrix of the random effects is replaced with a suitable proxy matrix. Here, the general idea used is as follows. After writing the joint density $f(y, \gamma)$ of the response variable $y$ and the random effects vector $\gamma$, they consider the penalized profiled likelihood $L_n(\beta, \hat{\gamma}(\beta)) - n \sum_{j=1}^{d_n} p_\lambda(|\beta_j|)$, where $L_n(\beta, \hat{\gamma}(\beta)) = f(y, \hat{\gamma}(\beta))$ with $\hat{\gamma}(\beta)$ the empirical Bayes estimate of $\gamma$ [88] in which the covariance matrix of $\gamma$ is replaced with the proxy matrix. $d_n$ may increase with the sample size $n$ and $p_\lambda$ is a concave penalty function (SCAD and LASSO in their simulation and application section). Although the proxy matrix may be different from the true one, it may still yield correct model selection results at the cost of some additional bias [54]. [164] deal with theoretical and computational aspects for high-dimensional selection of fixed effects in the linear mixed-effects model, where the consistency of the estimator is proven along with a non-asymptotic oracle result for the adaptive LASSO estimator, under the assumption that the eigenvalues of $Z^\top Z$ are bounded. Here, $Z$ is the random effects design matrix, and an explicit analytical expression of the regularization parameter $\lambda$ is also given.

Other approaches for variable selection in linear mixed-effects models include “fence” pro-
cedure and Bayesian techniques. [98] introduce a class of strategies known as a fence methods
intended for variable selection in both linear and generalized linear mixed-effects models. The
fence strategy is based on a measure of lack-of-fit that is a quantity $Q_M = Q_M(y, \theta_M)$, where $y$ is
the response variable, $M$ indicating a candidate model for the selection and $\theta_M$ denotes the vec-
tor of parameters under $M$. Counting among the typically rare and ad-hoc selection procedures,
the fence method is computationally very demanding, particularly because it involves the estimation of the standard deviation of the difference of lack-of-fit measures. For more details on these ad-hoc selection procedures, see the nice review paper of [131] and references therein. Bayesian model selection requires to assign a prior distribution over the model parameters and compute the posterior probabilities of each of them. These computations can be difficult so are usually carried out by applying sophisticated Markov Chain Monte Carlo (MCMC) algorithms [131]. Examples of Bayesian model selection include [33] and [160] who point out that these kinds of MCMC methods are generally time consuming to implement, requiring special software and depend on subjective choice of the hyperparameters in the priors.

In this paper, we discuss the selection of fixed effects in the linear mixed-effects model using an adaptive ridge (AR) penalty of the profiled log-likelihood, where the random effect covariance matrix is Cholesky factorized for solving a preliminary penalized least square problem. The profiled likelihood is calculated by slightly modifying the approach proposed by [14]. The weights matrix of the AR procedure introduced here is updated in such a way that the procedure converges toward selection with $L_0$ penalty, as have done [65]. The present selection strategy is intended to both low and high-dimensional settings.

The rest of the article is organized as follows. The profiled log-likelihood is introduced in Section 2. In Section 3, we present the weighted ridge procedure, and Section 4 presents the simulation studies.

## 4.2 profiled log-likelihood for the linear mixed-effects model

In this section, we consider the classical linear mixed-effect model setting where the number of observations $n$ is larger than the number of covariates $p$. By slightly modifying the approach introduced by [14], we calculate the profiled likelihood function.

### 4.2.1 Model and notations

We consider the linear mixed-effects model in which the residual terms are homoscedastic and independent of the random effects as follows.

$$ \mathcal{Y} = X\beta + Z\gamma + \varepsilon, \qquad (4.2) $$

$$ \gamma \sim \mathcal{N}(0, \Gamma), \quad \varepsilon \sim \mathcal{N}(0, \sigma^2 I_n) \quad \text{and} \quad \gamma \perp \varepsilon, \qquad (4.3) $$

where $\mathcal{Y}$ is the random response variable whose observed value is $y_{\text{obs}} \in \mathbb{R}^n$, $\gamma \in \mathbb{R}^q$ is the unobserved random effects vector with covariance matrix $\Gamma$, $\varepsilon$ is the residual term, $\beta \in \mathbb{R}^p$ is the fixed effects vector, and $X$ and $Z$ are the fixed and random effects related design matrices of dimensions $n \times p$ and $n \times q$, respectively. $\sigma^2 I_n$ is the covariance matrix of $\varepsilon$ with $\sigma > 0$ and $I_n$ the $n \times n$ identity matrix.

As a variance-covariance matrix, $\Gamma$ must be positive semidefinite. Conveniently, the model is expressed in terms of a relative covariance factor, $\Lambda_{\theta}$, which is a $q \times q$ matrix, depending on the variance components vector, $\theta$, that generate the symmetric $q \times q$ variance-covariance matrix, $\Gamma$,
according to

$$ \Gamma = \sigma^2 \Lambda_{\theta} \Lambda_{\theta}^{\top}, \qquad (4.4) $$

where $\sigma$ is the same scale parameter as in Equation (4.3). This factorization of $\Gamma$ yields the existence of a random vector $\mathcal{U}$ such that

$$ \gamma = \Lambda_{\theta} \mathcal{U}, \qquad (4.5) $$

with

$$ \mathcal{U} \sim \mathcal{N}(\mathbf{0}, \sigma^2 I_q) \qquad (4.6) $$

which is called a spherical random effects$^{(1)}$ variable. The model can therefore be re-expressed as follows

$$ \mathcal{Y} = X\beta + Z\Lambda_{\theta}\mathcal{U} + \varepsilon, \qquad (4.7) $$

$$ \mathcal{U} \sim \mathcal{N}(\mathbf{0}, \sigma^2 I_q), \quad \varepsilon \sim \mathcal{N}(\mathbf{0}, \sigma^2 I_n) \quad \text{and} \quad \mathcal{U} \perp \varepsilon. \qquad (4.8) $$

The parameters of the model are the fixed effects vector $\beta$ and the variance components $\theta$ and $\sigma^2$. This formulation of the linear mixed-effects model allows, not only, the use of a singular matrix $\Lambda_{\theta}$ which arises in practice, but also for a relatively compact expression for the profiled log-likelihood of $\beta$ and $\theta$, conditional on $y_{\text{obs}}$.

### 4.2.2 profiled likelihood

The profiled likelihood considered here is expressed through the following theorem.

**Theorem 4.2.1.** Suppose that $y$ is a realization of a random vector $\mathcal{Y}$ satisfying the linear mixed-effects model expressed by Equations (4.7) and (4.8), where $\beta$, $\theta$ and $\sigma^2$ are the parameters to be estimated. Denoting by $L_{\theta}$ the matrix such that $L_{\theta}^{\top} L_{\theta} = (Z\Lambda_{\theta})^{\top} Z\Lambda_{\theta} + I_q$, $\tilde{u}$ the conditional mean of $\mathcal{U}$ given that $\mathcal{Y} = y$, and $g(\tilde{u}) = \|y - X\beta - Z\Lambda_{\theta}\tilde{u}\|^2 + \|\tilde{u}\|^2$, the profiled log-likelihood of $\beta$ and $\theta$ conditional on $y$ is

$$ \tilde{\ell}(\beta, \theta|y) = -\frac{1}{2} \log |L_{\theta}|^2 - \frac{n}{2} \left[ 1 + \log \left( \frac{2\pi g(\tilde{u})}{n} \right) \right]. \qquad (4.9) $$

*Proof.* Denoting by $u$ a realization of $\mathcal{U}$, the density of $y$ is expressed as

$$ f(y) = \int_{\mathbb{R}^q} f(y, u)du = \int_{\mathbb{R}^q} f(y|u)f(u)du. $$

$$ f(y|u)f(u) = (2\pi\sigma^2)^{-(n+q)/2} \exp \left[ -\frac{\|y - X\beta - Z\Lambda_{\theta}u\|^2 + \|u\|^2}{2\sigma^2} \right], \qquad (4.10) $$

where $\|u\|^2 = u^{\top}u$, with $^{\top}$ denoting the transpose operator.

$$ \|y - X\beta - Z\Lambda_{\theta}u\|^2 + \|u\|^2 = \left\| \begin{pmatrix} y - X\beta \\ 0 \end{pmatrix} - \begin{pmatrix} Z\Lambda_{\theta} \\ I_q \end{pmatrix} u \right\|^2 = g(u), \qquad (4.11) $$

<footnote>$^{(1)}$[14] argued that the term “spherical” is related to the fact that the contours of the $\mathcal{N}(0, \sigma^2 I_q)$ probability density are spheres.</footnote>
and solving the penalized least squares problem that is to minimize $g(u)$ over $u$ implies

$$ \tilde{u} = \arg \min_{u \in \mathbb{R}^q} g(u) \iff \begin{pmatrix} Z\Lambda_{\theta} \\ I_q \end{pmatrix}^{\top} \begin{pmatrix} Z\Lambda_{\theta} \\ I_q \end{pmatrix} \tilde{u} = \begin{pmatrix} Z\Lambda_{\theta} \\ I_q \end{pmatrix}^{\top} \begin{pmatrix} y - X\beta \\ 0 \end{pmatrix}. \quad (4.12) $$

Viewed as a function of $u$, $g$ is $C^{\infty}$ and the first and the second differential of $g$ are

$$ \mathrm{d}g(u) = 2 \mathrm{tr} \left\{ -(y - X\beta - Z\Lambda_{\theta}u)^{\top} Z\Lambda_{\theta} \mathrm{d}u + u^{\top} \mathrm{d}u \right\} $$

and

$$ \mathrm{d}^2 g(u) = 2 \mathrm{tr} \left\{ (Z\Lambda_{\theta})^{\top} Z\Lambda_{\theta} \mathrm{d}u \mathrm{d}u^{\top} + \mathrm{d}u \mathrm{d}u^{\top} \right\}. $$

This implies that

$$ \frac{\partial g}{\partial u}(u) = -2 (y - X\beta - Z\Lambda_{\theta}u)^{\top} Z\Lambda_{\theta} + 2u^{\top} $$

and

$$ \frac{\partial^2 g}{\partial u \partial u^{\top}}(u) = 2 (Z\Lambda_{\theta})^{\top} Z\Lambda_{\theta} + 2I_q. $$

Then, $g(u)$ can be rewritten as

$$ g(u) = g(\tilde{u}) + \frac{\partial g}{\partial u}(\tilde{u})(u - \tilde{u}) + \frac{1}{2}(u - \tilde{u})^{\top} \frac{\partial^2 g}{\partial u \partial u^{\top}}(\tilde{u})(u - \tilde{u}). \quad (4.13) $$

Referring to Equation (4.12),

$$ \begin{pmatrix} Z\Lambda_{\theta} \\ I_q \end{pmatrix} \left[ \begin{pmatrix} y - X\beta \\ 0 \end{pmatrix} - \begin{pmatrix} Z\Lambda_{\theta} \\ I_q \end{pmatrix} \tilde{u} \right] = 0 \implies \frac{\partial g}{\partial u}(\tilde{u}) = 0. \quad (4.14) $$

Then,

$$ \begin{align} g(u) &= g(\tilde{u}) + \frac{1}{2}(u - \tilde{u})^{\top} \frac{\partial^2 g}{\partial u \partial u^{\top}}(\tilde{u})(u - \tilde{u}) \nonumber \\ &= g(\tilde{u}) + (u - \tilde{u})^{\top} [(Z\Lambda_{\theta})^{\top} Z\Lambda_{\theta} + I_q] (u - \tilde{u}) \nonumber \\ &= g(\tilde{u}) + \|L_{\theta}(u - \tilde{u})\|^2, \quad \text{with } L_{\theta}^{\top} L_{\theta} = (Z\Lambda_{\theta})^{\top} Z\Lambda_{\theta} + I_q, \tag{4.15} \end{align} $$

and

$$ \begin{align} f(y) &= (2\pi\sigma^2)^{-(n+q)/2} \int_{\mathbb{R}^q} \exp \left[ -\frac{g(\tilde{u}) + \|L_{\theta}(u - \tilde{u})\|^2}{2\sigma^2} \right] du \nonumber \\ &= (2\pi\sigma^2)^{-(n+q)/2} \exp \left[ -\frac{g(\tilde{u})}{2\sigma^2} \right] \int_{\mathbb{R}^q} \exp \left[ -\frac{\|L_{\theta}(u - \tilde{u})\|^2}{2\sigma^2} \right] du. \tag{4.16} \end{align} $$

$$ v = L_{\theta}(u - \tilde{u}) \implies du = \frac{1}{|L_{\theta}|} dv, \quad \text{and} $$
$$ 
\begin{align} 
f(y) &= (2\pi\sigma^2)^{-n/2} |L_{\theta}|^{-1} \exp\left[-\frac{g(\tilde{u})}{2\sigma^2}\right] \int_{\mathbb{R}^q} (2\pi\sigma^2)^{-q/2} \exp\left[-\frac{\|v\|^2}{2\sigma^2}\right] dv \nonumber \\ 
&= (2\pi\sigma^2)^{-n/2} |L_{\theta}|^{-1} \exp\left[-\frac{g(\tilde{u})}{2\sigma^2}\right] \tag{4.17} 
\end{align} 
$$

The log-likelihood of $\beta$, $\theta$ and $\sigma^2$ conditional on $y$ is

$$ \ell(\beta, \theta, \sigma^2|y) = -\frac{n}{2}\log(2\pi\sigma^2) - \frac{1}{2}\log|L_{\theta}|^2 - \frac{g(\tilde{u})}{2\sigma^2}. \quad (4.18) $$

$\frac{\partial \ell(\beta, \theta, \sigma^2|y)}{\partial \sigma^2} = 0 \implies \sigma^2 = \frac{g(\tilde{u})}{n}$. By profiling out $\sigma^2$, the profiled log-likelihood, $\tilde{\ell}(\beta, \theta|y)$, of $\beta$ and $\theta$ conditional on $y$ is expressed through

$$ -2\tilde{\ell}(\beta, \theta|y) = \log|L_{\theta}|^2 + n \left[ 1 + \log \left( \frac{2\pi g(\tilde{u})}{n} \right) \right]. $$

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

## 4.3 $L_0$ estimator of $\beta$ using iteratively weighted ridge procedure

The $L_0$ estimator of the fixed-effects vector $\beta$ using iteratively weighted ridge procedure presented here fits both the low and the high-dimensional settings. Theoretically, this selection method may enjoy great stability since it uses neither inverse of $X$ (fixed-effects design matrix) nor inverse of $Z$ (random effects design matrix).

### 4.3.1 Adaptive Ridge penalty for the profiled likelihood

Let us assume that the true underlying fixed-effects vector $\beta_{\text{true}}$ is sparse in the sense that many of its coefficients are zero. To enforce the sparsity of the estimator of $\beta$, we add a weighted $L_2$ (ridge) penalty for the fixed-effects vector $\beta$ to the profiled log-likelihood function. Thus, we are considering the objective function

$$ \tilde{\ell}_{\lambda,w}(\beta, \theta) = -2\tilde{\ell}(\beta, \theta|y_{\text{obs}}) + \lambda\beta^\top W\beta, \quad (4.19) $$

where $W = \text{diag}(w_1, \dots, w_p)$ is a $p \times p$ diagonal matrix and $\lambda \ge 0$ is a regularization parameter. We aim at estimating $\beta$, $\theta$ and $\sigma^2$ by

$$ (\tilde{\beta}_{\lambda,w}, \tilde{\theta}_{\lambda,w}) = \arg\min_{\beta, \theta} \tilde{\ell}_{\lambda,w}(\beta, \theta) \quad \text{and} \quad \tilde{\sigma}^2_{\lambda,w} = \frac{\|y - X\tilde{\beta}_{\lambda,w} - Z\Lambda_{\tilde{\theta}_{\lambda,w}}\tilde{u}\|^2 + \|\tilde{u}\|^2}{n}, \quad (4.20) $$

where $\tilde{u}$ is calculated by the penalized least squares algorithm that has been hinted at Equation (4.12). The criterion $\tilde{\ell}_{\lambda,w}(\beta, \theta)$ is minimized using one of the constrained optimization functions in $\mathbb{R}$ [146], to provide the estimators $\tilde{\beta}_{\lambda,w}$, $\tilde{\theta}_{\lambda,w}$ and $\tilde{\sigma}^2_{\lambda,w}$.
### 4.3.2 Iteratively Weighted ridge procedure

The $L_0$ penalty for regularization arises naturally in many classical model selection since, indeed, it counts the number of non-vanishing parameters, giving a nice interpretation of the best subset selection and admits nice sampling properties [11]. However, its computation is infeasible in high dimensional settings and clearly argued to be a combinational problem with NP-complexity [53]. Some workarounds are reported in the literature, where the proposed procedure converge toward the $L_0$-penalty based selection. For example, [65] introduced an adaptive ridge procedure that helps to approximate $L_0$-penalty performances. Here, we are using the same procedure where the weight matrix ($W$ defined in Equation (4.19)) diagonal elements $w_j$, $1 \le j \le p$ are iteratively computed and defined as

$$w_j^{(k)} = \left[ |\beta_j^{(k)}|^2 + \delta^2 \right]^{-1}, \qquad (4.21)$$

as have done [65], taking inspiration from [20, 28, 78, 154] and their simulation results. In Equation (4.21), $k$ identifies the iteration and $\beta_j$ is the $j$th component of $\beta$. The general expression of $w_j^{(k)}$ is $w_j^{(k)} = \left[ |\beta_j^{(k)}|^\tau + \delta^\tau \right]^{\frac{q-2}{\tau}}$, where $q$ precises the norm $\| \cdot \|_{L_q}$ for the penalty, $\delta$ calibrates which effect sizes are considered relevant and $\tau$ determines the quality of the approximation $w_j \beta_j^2 \approx |\beta_j|^q$. In practice, $\delta = 10^{-5}$ seems to perform well. For more details, see [65]. More precisely, the selection procedure performed here is as follows. For a fixed $\lambda$, $\tilde{\beta}_{\lambda,w}$ is initialized at $(1, \dots, 1)^\top$, the vector selection, say, which identifies the selected fixed effects is initialized at $(1, \dots, 1)^\top$ and $W$ is initialized at $\text{diag}(1, \dots, 1)$. Then the steps come.

1) $\tilde{\beta}_{\lambda,\text{old}} \leftarrow \tilde{\beta}_{\lambda,w}$

2) perform the optimization $(\tilde{\beta}_{\lambda,w}, \tilde{\theta}_{\lambda,w}) = \arg \min_{\beta, \theta} \tilde{\ell}_{\lambda,w}(\beta, \theta)$, initializing $\beta$ by $\tilde{\beta}_{\lambda,\text{old}}$ and $\theta$ by $\theta_0$. Here, the components of $\theta$ which are variances are initialized by 1 and those which are not variances are initialized by 0. Thus, $\theta_0$ components are 0 or 1.

3) $w_j \leftarrow (\tilde{\beta}_{\lambda,w,j}^2 + \delta^2)^{-1}$, where $\tilde{\beta}_{\lambda,w,j}$ is the $j$th component of $\tilde{\beta}_{\lambda,w}$, and $w_j$ is the $j$th element of $W$'s diagonal.

4) selection$_{old} \leftarrow selection$ and selection $\leftarrow W \cdot \text{diag}(\tilde{\beta}_{\lambda,w}) \cdot \tilde{\beta}_{\lambda,w}$

5) if $|\text{selection} - \text{selection}_{\text{old}}| < \text{tol} = 10^{-5}$, then the selection can be considered as well performed for $\lambda$. Thus, we choose a new value for $\lambda$. If $|\text{selection} - \text{selection}_{\text{old}}| \ge \text{tol}$, the selection does not perform well and we go to the item 1) by choosing the current $\tilde{\beta}_{\lambda,w}$ as $\tilde{\beta}_{\lambda,\text{old}}$, without changing the value of $\lambda$.

For some $\lambda$ values, the selection vector may contain 0 or 1 as components. If it is 1, the corresponding fixed effect $\beta_j$ is selected, and if it is 0 the corresponding $\beta_j$ is not selected for $\lambda$.

For the choice of the regularization parameter $\lambda$, we propose to use the Bayesian Information Criterion (BIC) criterion defined as

$$c_{n,\lambda} = -2\ell(\hat{\beta}_{\lambda,\text{sel}}, \hat{\theta}_{\lambda,\text{sel}}, \widehat{\sigma}_{\lambda,\text{sel}}^2|y_{\text{obs}}) + \log(n) \cdot \hat{d}_{\lambda}, \qquad (4.22)$$
where $\hat{\beta}_{\lambda,sel}$, $\hat{\theta}_{\lambda,sel}$ and $\hat{\sigma}^2_{\lambda,sel}$ are the ML parameters estimators considering the selected variables. $\hat{d}_{\lambda} = \#\{\hat{\beta}_{\lambda,sel,j} \neq 0 : 1 \le j \le p\} + \dim(\theta) + 1$ is the sum of the number of nonzero fixed-effects and the number of variance components. This form of $\hat{d}_{\lambda}$ has been suggested by [17] and empirically validated by [165].

## 4.4 Simulation studies

In this section, we assess the performance of our approach using simulated data. We compare the obtained results with those coming from the Lasso implementation in a situation where the number of noise variables are excessive.

Since we use one of the optimizers available in the R software for minimizing the $\tilde{\ell}_{\lambda,w}(\beta, \theta)$ criterion, for too bigger values of $p$, especially when $n \ll p$, the convergence of the used algorithm (“nlminb” for example) is hardly or no more reached. Due to this convergence problem, we restrict the simulation studies to the low-dimensional setting where $p > 40$ and will focus on this problem in another coming paper with the same theoretical approach.

We restrict ourselves to the case of longitudinal data study with $N$ observations coming from $n$ subjects where each subject $i$ has $n_i$ observations. The “working” data sets are simulated under the following model.

$$y_i = X_{1i}\beta^* + \gamma_i 1_{n_i} + \varepsilon_i, \quad \gamma_i \sim \mathcal{N}(0, \Gamma), \quad \varepsilon_i \sim \mathcal{N}(0, \sigma^2 I_{n_i}), \quad \text{for } i = 1, \dots, n; \quad (4.23)$$

where $\beta^* \in \mathbb{R}^{p_1}$ is the true fixed effects vector, $X_1$ is a $N \times p_1$ design covariates matrix, $1_{n_i} = (1, \dots, 1) \in \mathbb{R}^{n_i}$. $\gamma \perp \varepsilon$ and $\gamma_i \perp \gamma_{i'}$ for $i \neq i'$. $\gamma_i$ is a random intercept for the $i$th subject and $\varepsilon$ is the residual term of the model.

We suppose that we are following up a sample of subjects where the goal is to evaluate how their weights are influenced by other variables including the age, the sex and the nutrition score “nscore”, and which variables govern this influence. The covariates sex, nscore and age are staked in the model matrix $X_1$ and all other covariates are staked in another $N \times p_2$ model matrix $X_2$ such that $X = (X_1|X_2)$ with $\dim(X) = N \times p$ and $p = p_1 + p_2$. The components of the vector (variable) age are randomly sampled from a uniform distribution in [18, 37] and the nscore variable is also uniform in [20, 50]. We would like to fit to the data, the model

$$y_i = X_i\beta + \gamma_i 1_{n_i} + \varepsilon_i, \quad \gamma_i \sim \mathcal{N}(0, \Gamma), \quad \varepsilon_i \sim \mathcal{N}(0, \sigma^2 I_{n_i}), \quad \text{for } i = 1, \dots, n; \quad (4.24)$$

where $\beta = (\beta_1^\top, \beta_2^\top)^\top$, with $\beta_1 \in \mathbb{R}^{p_1}$ and $\beta_2 = 0 \in \mathbb{R}^{p_2}$. We are then challenging to identify which $\beta$ components are zero.

For the working data sets, we choose $p_1 = 4$, $p = 54$, $N = 300$, $n = 90$, $\sigma = 1$, $\Gamma = 1$ and $\beta^* = (1, -1, -1, 1)$. For the fixed-effects $\beta$, we have in fact fifty zeros components and only four components are not zeros. We simulate 100 data sets for which we perform both the lasso selection and the IWR (Iteratively Weighted Ridge) procedure. One hundred values of the regularization parameter $\lambda$ are chosen in $[10^{-2}, 10^2]$. For each replication, the $\lambda$ value that minimizes the BIC criterion is retained for selecting the significant fixed effects. Denoting by $\beta^{**} = (1, -1, -1, 1, 0, \dots, 0)$, the mean squared error MSE = $E[\|\hat{\beta} - \beta^{**}\|_2^2]$ (with $\hat{\beta}$ coming
from lasso or IWR) is computed over the 100 replications. The cardinality of the estimated active set (i.e. $|S(\hat{\beta})|$, with $S(\hat{\beta}) = \{k : \hat{\beta}_k \neq 0\}$) is computed as well as the proportion TP of true positive (i.e. the selected set is exactly the true one). We also compute the proportion TPC in which the selected set contains the true one and the proportion ZP in which the true zeros are estimated.

<table><caption>Table 4.1: Selection performances comparison between LASSO and IWR procedures.</caption><thead><tr><th>Performance criterion</th><th>LASSO</th><th>IWR</th></tr></thead><tbody><tr><td>MSE</td><td>1.200</td><td>0.254</td></tr><tr><td>|S(β̂)|</td><td>4(1.614)</td><td>5(1.250)</td></tr><tr><td>TP</td><td>0%</td><td>35%</td></tr><tr><td>TPC</td><td>16%</td><td>90%</td></tr><tr><td>ZP</td><td>98%(0.028)</td><td>98%(0.023)</td></tr></tbody></table>

The selection results based on the 100 simulated data sets are indeed summarized through five performance criteria that are contained in Table 4.1, where the numbers between parentheses are the standard deviations related to the criteria mean values (printed just before these parentheses). Information from Table 4.1 show that the IWR procedure outperforms the LASSO one. IWR selects the correct model 35% of the time when LASSO has never found it (see TP values in Table 4.1). For information, we use the R software package lmllasso [163] to perform LASSO selection. 90% of the time, the model selected by IWR contains the true one against 16% for LASSO (see TPC values from Table 4.1). Obviously, the estimations are of better qualities from IWR than from LASSO (MSE = 0.254 for IWR, and MSE = 1.200 for LASSO).

The true zero estimation proportions (ZP) are computed in Table 4.2 which contains also their number of occurrence over the simulated data sets. The LASSO seems to find more often than

<table><caption>Table 4.2: True zero estimation proportions (ZP) with the number of occurrence over the 100 replications.</caption><thead><tr><th></th><th colspan="7">LASSO</th><th colspan="6">IWR</th></tr></thead><tbody><tr><td>ZP</td><td>0.88</td><td>0.90</td><td>0.92</td><td>0.94</td><td>0.96</td><td>0.98</td><td>1</td><td>0.90</td><td>0.92</td><td>0.94</td><td>0.96</td><td>0.98</td><td>1</td></tr><tr><td>Number</td><td>2</td><td>2</td><td>3</td><td>6</td><td>14</td><td>24</td><td>49</td><td>2</td><td>2</td><td>8</td><td>16</td><td>32</td><td>40</td></tr></tbody></table>

IWR all the true zeros (49% for LASSO against 40% for IWR, in Table 4.2). This may explain the overfitting behavior of IWR ($|S(\hat{\beta})| = 5$ for IWR and $|S(\hat{\beta})| = 4$ for LASSO in Table 4.1). It therefore seems that LASSO has a stronger shrinkage capability than has IWS which shows in turn a somewhat overfitting behavior than LASSO. Through these simulations studies, it appears that the iteratively weighted aspect of the IWR procedure higher the shrinkage performance of the ordinary Ridge and results in a selection method having better performance than the LASSO.
Like in the case of LASSO, it is possible for IWR to take advantage of a warm start of the
algorithm to obtain the full regularization path of the selection problem. For instance, Figure 4.1

<figure><img src="image_12.png" /><figcaption>Figure 4.1: Example of a full regularization path for the iteratively weighted ridge selection procedure with $N = 300$, $n = 90$, $p = 54$, $\sigma = 1 = \Gamma$. Only the first four components of $\beta^*$ are nonzero. The dataset is simulated using the R software with set.seed(3). The vertical red dashed bar corresponds to the minimum of BIC and shows the selected variables.</figcaption></figure>

shows the full regularization path from one of the simulated data. The vertical red dashed bar
corresponds to the minimum of BIC and shows the selected variables. On Figure 4.1, we clearly
have $|S(\hat{\beta})| = 4$, i.e., four variables selected at the end of the procedure.

**4.5 Conclusion**

In this paper, we have focused on the fixed-effects selection problem in the linear mixed-effects
model. We have introduced an iteratively weighted ridge procedure which enhances the shrinkage
performance of the ordinary ridge in order to approximate the performance of the $L_0$ based penalty
selection. This selection method is based on an adaptive ridge penalty of the profiled likelihood,
where the covariance matrix of the random effects is Cholesky factorized. The procedure fits both
the low and the high-dimensional settings and may enjoy a great numeric stability since it needs
no use of the inverse of the fixed-effects design matrix or the design matrix of the random-effects.
Due the problems of the lack of convergence in the high-dimensional settings using available
optimizers, we restrict the simulations studies to the low-dimensional case where we find out
that our selection procedure outperforms the LASSO selection. In another forthcoming paper, we
will focus on this convergence problem in the high-dimensional cases with the same theoretical
approach.
# APPLICATION À L'ÉTUDE DE L'ACQUISITION IMMUNITAIRE CONTRE LE PALUDISME CHEZ L'ENFANT À TORI-BOSSITO (BÉNIN)

## 5.1 Paludisme

Le paludisme est une maladie infectieuse à transmission vectorielle faisant intervenir trois acteurs que sont l'Homme, jouant le rôle d'hôte, infecté par un protozoaire (le parasite) du genre *Plasmodium* qui lui a été transmis par la piqûre d'un moustique (le vecteur) femelle du genre *Anopheles*. Outre les facteurs biologiques de l'hôte et du parasite, les facteurs géographiques sont essentiels pour comprendre la transmission du paludisme et expliquer la susceptibilité des individus aux infections. Le nombre de piqûres reçues par un individu dépend non seulement de l'abondance de gîtes larvaires dans l'environnement, mais aussi de la proportion d'anophèles vecteurs et d'anophèles porteurs de l'agent pathogène [43]. L'infectiosité des anophèles dépend directement du nombre de parasites circulants et par conséquent de la durée du cycle parasitaire qui elle-même dépend de la température du milieu de vie du vecteur. Par exemple, pour le *P. falciparum*, en dessus de l'optimum de 25°C de température, le cycle s'allonge et par conséquent la transmission diminue. Ainsi, la température, la pluviométrie et la végétation sont des facteurs essentiels à l'explication de la transmission palustre [37].

### 5.1.1 Vecteurs

Les vecteurs du paludisme humain appartiennent tous au genre *Anopheles* qui est de la famille des *Culicidae*, de l'ordre des *Dieptra*. Il existe plus de 450 espèces d'anophèles recensées, regroupées par régions géographiques et dont seulement 70 à 80 sont considérées comme étant des vecteurs. En Afrique tropicale, les cinq principaux vecteurs sont : *An. gambiae s.s.*, *An. aarabiensis*, *An. funestus*, *An. nili s.l.* et *An. moucheti*. En milieu tropical, un anophèle mâle a une durée de vie d'environ 10 jours alors que la femelle quant à elle peut vivre jusqu'à deux à quatre semaines. La femelle fécondée se nourrit de jus sucré et de sang qu'elle prélève tous les deux à trois jours sur un hôte vertébré. Elle pond séparément, à la surface de l'eau, 40 à 100 œufs qui éclosent au bout de 24 à 48 heures selon la température, et part ensuite à la recherche d'un nouvel hôte pour un autre repas sanguin [130].
5.1.2 Agent pathogène

Lors de son repas sanguin, l'anophèle libère dans le sang humain les parasites responsables du paludisme chez l'Homme. Ces parasites sont des protozoaires du genre *Plasmodium*. Cinq espèces de *Plasmodium* peuvent parasiter l'Homme : *P. malariae*, *P. vivax*, *P. knowlesi*, *P. ovale* et *P. falciparum* qui est le plus mortel, présentant une large gamme de manifestations pathologiques. Dans les régions équatoriales, il est transmis toute l'année avec des recrudescences saisonnières et ne survient qu'en période chaude et humide dans les régions sub-tropicales. Environ trois-quarts des cas d'infection due au *P. falciparum* dans le monde sont enregistrés en Afrique [180]. Ces dernières années, il a été enregistré certains cas humains de paludisme à *P. knowlesi*, un paludisme du singe rencontré dans certaines zones de forêts d'Asie du Sud-Est.

Le cycle biologique du *Plasmodium* se décompose en une phase sexuée chez l'Anopheles (hôte définitif) et en une phase asexuée chez l'Homme (hôte intermédiaire). Les *Plasmodium* pénètrent dans l'organisme humain sous forme de sporozoites à la faveur d'une piqûre d'un moustique infec-tieux. Les sporozoites se multiplient ensuite dans les cellules du foie après y avoir été transportés par la circulation sanguine (Voir Figure 5.1). Libérés dans le sang sous forme de mérozoites, ils envahissent les globules rouges. Il en résulte des accès de fièvre. A maturité, le schizonte hépatique éclate, libérant des mérozoites, formes uninucléées qui poursuivront leur développement au cours de la phase érythrocytaire. L'invasion de l'érythrocyte débute par la liaison du mérozoite à la surface d'un érythrocyte. Les schizontes peuvent ainsi infecter d'autres globules rouges ou se transformer en gamétocytes mâles et femelles. Par la suite, un moustique se contamine par piqûre en absorbant du sang contenant des gamétocytes qui passeront à l'étape de gamètes dans le tube digestif de l'insecte. La fécondation d'un gamète femelle par un mâle produit un zygote (cel-

<figure><img src="image_13.png" /><figcaption>FIGURE 5.1 : Cycle du développement et de la reproduction des *Plasmodium*. (CDC, 2007)</figcaption></figure>

lule œuf) qui se développe en sporozoite. Les sporozoites migrent ensuite dans les glandes sali-
vaires du moustique, d'où ils pourront contaminer un nouvel individu lors d'une piqûre. Ainsi, les éléments asexués (trophozoites et schizontes) sont digérés après absorption des différents stades du parasite par le moustique pendant le repas sanguin ; seuls les gamétocytes poursuivront leur développement. Rapidement, le gamète femelle se transforme en macrogamète et le gamétocyte mâle donne naissance à 8 microgamètes flagellés qui se dirigent à la rencontre du macrogamète. La fécondation donne naissance à un œuf mobile (l'ookinète) qui traverse la paroi de l'estomac, formant l'oocyte dans lequel s'individualisent les sporozoites qui, à leur tour, gagneront les glandes salivaires de l'anophèle [41].

### 5.1.3 Accès palustre grave et groupes à risques

*P.falciparum* est a l'origine des accès palustres sévères qui sévissent essentiellement chez les individus non immuns et certains groupes à risques parmi les populations vivant en zone d'endémie - femmes enceintes et enfants de moins de cinq ans - [43]. Chez un malade présentant une parasitémie avec des formes asexuées de *P.falciparum* sans aucune cause manifeste de symptômes, la présence de troubles de la conscience, de prostration, de convulsions multiples, d'ictère clinique accompagné d'autres signes de dysfonctionnement des organes vitaux,... indique qu'il souffre d'un paludisme grave [166].

Les femmes enceintes, particulièrement pendant leur première grossesse, sont susceptibles aux nouvelles infections palustres et montrent une séquestration massive d'hématies parasitées dans l'espace sanguin maternel du placenta. Or le paludisme gestationnel est généralement asymptomatique mais peut avoir de sérieuses conséquences sur le devenir de la mère et du fœtus : anémie maternelle, faible poids de naissance, mort fœtale [43, 69]. Les enfants de plus de cinq ans qui

<figure><img src="image_14.png" /><figcaption>FIGURE 5.2 : Mère et enfant face au paludisme (IRD, 2010).</figcaption></figure>

ont été exposés régulièrement à P. falciparum développent une immunité dite "non-stérilisante" ou
immunité “clinique” qui n’apparaît qu’après un contact continu avec une large variété d’antigènes
parasitaires [100]. Cette semi-immunité permet au patient d’être protégé contre la déclaration de
symptômes, mais ne permet pas d’éliminer les parasites infectants [43, 199].

5.1.4 A la recherche d'un vaccin contre le paludisme

Le *P. falciparum* est connu pour son polymorphisme et sa variabilité génétique qui est relative aux protéines exprimées à la surface des érythrocytes - variant de surface antigénique (VAR) - codés par des gènes ayant plusieurs copies. Le polymorphisme quant à lui concerne les gènes représentés en une seule copie dans le génome et dépend de l'étendue des allèles disponibles et de la reproduction sexuée. D'après les études menées sur les gènes *var* - l'un des gènes les plus polymorphes - et le phénotype de cytoadhésion du parasite, le taux de mutation estimé varie entre 2 et 18 par génération [71]. Lors d'une infection multi-clonale, la majeur partie des anticorps sont synthétisés spécifiquement contre le VAR le plus représenté. Les parasites porteurs de cet antigène seront neutraлизated alors que les parasites mutants seront épargnés, pouvant ainsi se développer et amener le système immunitaire à produire à nouveau des anticorps spécifiques du VAR majoritaire. Et ces parasites échappent au système immunitaire. Cette variabilité est donc à la base de toute la difficulté du système immunitaire à éliminer le parasite, et par conséquent du réel défi que représente la mise en place d'un vaccin efficace sur une longue durée [43]. Par ailleurs, le polymorphisme allélique de *P. falciparum* explique pourquoi il est nécessaire que le l'Homme soit exposé de nombreuses fois à divers antigènes avant d'acquérir une semi-immunité prévenant les accès palustres symptomatiques [100].

Trois stades dans le cycle de *P. falciparum* peuvent être la cible d'un vaccin : le stade pré-érythrocytaire, le stade érythrocytaire et le stade sexué [43, 77]. Les essais vaccinaux basés sur le stade pré-érythrocytaire ont pour objectif d'empêcher toute libération de mérozoites provenant du schizontes dans le sang. Des exemples de vaccins utilisant cette approche ont montré leur limite en phase IIb d'essai vaccinal, puisque incapables d'induire une immunité stérilisante [5, 179].

Les vaccins visant à stopper le développement des parasites asexués au stade érythrocytaire cherchent à inhiber l'invasion des hématies dans le but de contrôler la parasitémie et éviter ainsi l'évolution de l'infection vers des formes cliniques [43]. Ici, les candidats vaccins sont les antigènes portés par les mérozoites libérés par les schizontes érythrocytaires, essentiellement les Merozoite Surface Protein (MSP)-1, MSP2, MSP3, Apical Membrane Antigen (AMA)-1, Glutamate Rich Protein (GLURP). Des expériences de transfert d'immunoglobulines (Ig) purifiées d'adultes hyper-immuns à des enfants ont montré que les réponses anticorps jouaient un rôle dans l'acquisition de l'immunité adaptative dirigée contre les stades sanguins de *P. falciparum* [27, 34, 158].

Les vaccins dirigés contre le stade sexué de *P. falciparum* permettraient de bloquer la transmission vectorielle de la maladie, où les antigènes des gamètes en sont la cible. Ceci dit, il serait nécessaire de vacciner simultanément toute la population afin de stopper la transmission. De plus, cette vaccination devrait potentiellement être combinée à un vaccin des stades pré-érythrocytaire ou érythrocytaire afin d'espérer son efficacité.

La complexité du parasite, sa grande diversité antigénique et ses stratégies d'échappement immunitaire rendent le développement d'un vaccin difficile. Des efforts sont actuellement menés pour mettre au point un vaccin spécifiquement destiné à protéger la femme enceinte et son fœ-
tus aussi bien que le jeune enfant [43]. Les données que nous utilisons dans cette thèse sont is-
sues des recherches menées en vue d'étudier la réponse anticorps spécifique d'antigènes du stade
érythrocytaire, stade à l'origine des manifestations cliniques du paludisme.

## 5.2 Données pour l'acquisition de la réponse anticorps spécifique du paludisme

Les données que nous analysons dans le cadre de cette thèse ont été recueillies dans la commune de Tori-Bossito située dans le Nord-Ouest de Cotonou (au Bénin) à 40 km du centre ville (voir Figure 5.3). Cette commune, comptant quarante-sept villages et quartiers de ville, est caractérisée

<figure><img src="image_15.png" /><figcaption>FIGURE 5.3 : Localisation géographique de Tori-Bossito et les neuf villages d'étude au Bénin [108].</figcaption></figure>

par un climat sub-équatorial avec deux saisons de pluie et deux saisons sèches. La transmission palustre y est relativement identique chaque année (paludisme stable) fortement influencée par les pluies. Les vecteurs du paludisme rencontrés le plus souvent dans cette zone sont les *Anopheles gambiae* et les *Anopheles funestus* [35, 48].

Pour des raisons d'efficacité dans le suivi des sujets prenant part à l'étude, neuf villages parmi ceux qui sont proches des trois centres de santé de la commune de Tori-Bossito ont été sélectionnés. Il s'agit de : Avamé, Anavié, Cada-centre, Dohinoko, Gbédjougo, Gbétaga, Houngo, Zébé et Zoungoudo. Le programme avait pour objectif d'étudier les déterminants des premières infections palustres chez des nouveau-nés au cours des 18 premiers mois de vie. En outre il s'agissait également d'étudier la mise en place de la réponse immune spécifique des formes asexuées de *P. falciparum*. Pour atteindre ces deux objectifs il a été mis en place une cohorte de femmes
enceintes recrutées à l'accouchement et de leurs enfants suivis jusqu'à 18 mois. Etait susceptible d'être prise en compte par l'enquête, toute femme enceinte habitant l'un des villages retenus et accouchant dans l'une des trois maternités de la commune entre juin 2007 et juillet 2008. 620 femmes au total ont été incluses et chaque nouveau-né était suivi pendant ses dix-huit premiers mois de vie.

Les informations recueillies incluent des attributs descripteurs de la mère, de son nouveau-né aussi bien que de leur milieu de vie. A l'inclusion, était recueillies des informations épidémiologiques, comportementales et familiales, et celles nutritionnelles, environnementales, relatives à la malaria étaient enregistrées tout au long du suivi - voir [108] pour plus de détails. L'objectif de cette étude était d'évaluer les déterminants de l'incidence palustre dans les premiers mois de vie de l'enfant. Entre autres objectifs spécifiques, il s'agissait de comprendre comment est-ce que le système immunitaire de l'enfant va construire la capacité à répondre à des antigènes du stage érythrocytaire. Ainsi, les sept antigènes sélectionnés : AMA1, MSP1, MSP2 (3D7, FC27), MSP3 et GLURP (R0, R2) font partie de la phase érythrocytaire, sont immunogènes et sont des candidats vaccins prometteurs induisant une réponse anticorps protectrice démontrée [43]. Par ailleurs, ils sont localisés à la surface du mérozoite et sont donc accessibles aux anticorps lorsque les mérozoites sont libérés dans le sang. De plus, des protéines recombinantes les reproduisant sont disponibles. Pour plus de détails spécifiques à ces antigènes, voir par exemple [43].

**5.3 Applications**

Nous exposons ici deux illustrations relatives à l'analyse des données d'anticorps palustres mais que nous gardons sous anonymat pour des raisons de confidentialité. Ces données proviennent de l'UMR 216 et ont été codées et nous n'indiquons en conséquence aucun résultat biologique issus de nos analyses. Deux protéines recombiantes du *P. falciparum* ont été utilisées pour obtenir deux sous classes de IgG (IgG1 et IgG3) donnant ainsi des anticorps que nous avons nommés

$$
IgG1\_A1, IgG3\_A1, IgG1\_A2, IgG3\_A2, IgG1\_B, IgG3\_B, IgG1\_C, IgG3\_C, \quad (5.1)
$$

où A1, A2, B et C sont des antigènes du stade érythrocytaire. L'objectif de l'analyse de ces données dans le cadre de nos travaux de thèse est d'évaluer l'effet de l'infection palustre sur l'acquisition immunitaire de l'enfant pendant ces dix-huit premiers mois de vie. Puisque les antigènes qui caractérisent le statut immun de l'enfant inter-agissent entre eux dans l'organisme de l'enfant, nous analysons les caractéristiques de la distribution jointe de ces anticorps conditionnellement à l'infection palustre et à d'autres facteurs d'intérêt. Dans les modèles que nous ajustons à ces données, les variables dépendantes conc.Y (voir tableau 5.1) décrivent le niveau d'expression de l'antigène Y chez l'enfant à 3, 6, 9, 12, 15 et 18 mois, où Y désigne une protéine de la liste fournie à l'équation (5.1). Par ailleurs, toutes les autres variables (à l'exception de conc.Y) contenues dans le tableau 5.1 sont des covariables. Nous sommes ainsi en présence de huit variables dépendantes qui décrivent le profile longitudinale (chez l'enfant) des protéines listées à l'équation (5.1).
<table><caption>TABLE 5.1 : Variables prises en compte par nos analyses</caption><thead><tr><td>Variable</td><td>Description</td></tr></thead><tbody><tr><td>id</td><td>Identifiant de l'enfant</td></tr><tr><td>conc.Y</td><td>concentration de Y</td></tr><tr><td>conc_CO.Y</td><td>Concentration de Y mesurée dans le sang du cordon</td></tr><tr><td>conc_M3.Y</td><td>Concentration de Y prédite dans le sang périphérique de l'enfant à 3 mois</td></tr><tr><td>ap</td><td>Apposition placentaire</td></tr><tr><td>hb</td><td>Niveau d'hémoglobine</td></tr><tr><td>inf_trim</td><td>Nombre d'infections palustres dans les 3 derniers mois</td></tr><tr><td>pred_trim</td><td>Nombre moyen des prédictions entre 3-6 mois en matière d'exposition au moustique</td></tr><tr><td>nutri_trim</td><td>Moyenne des valeurs de nutrition mensuelle du trimestre écoulé</td></tr></tbody></table>

### 5.3.1 Première illustration : Classification hiérarchique de protéines palustres

Dans les modèles que nous avons ajustés aux données, nous avons mis un intercept aléatoire par enfant et une pente aléatoire par enfant suivant l'infection palustre (i.e. la variable inf_trim). L'illus-

<figure><img src="image_16.png" /><figcaption>FIGURE 5.4 : Arbre de classification hiérarchique de protéines palustres</figcaption></figure>

tration que nous faisons ici est d'analyser conjointement chacune des vingt-huit paires de protéines en vue de vérifier si certains profiles de protéines peuvent être analysés de façon indépendante, conditionnellement à la configuration du modèle ajusté. Après avoir mis en œuvre sur toutes ces paires de protéines le test de corrélation bi-varié introduit au Chapitre 2, en présence d'une correc- tion de Bonferroni, les p-valeur obtenues vont de 0.000 à 0.932. La p-valeur 0.932 est d'ailleurs la seule qui soit en dessous du seuil de 5% et provient de la paire de protéines (IgG3_A1, IgG1_B). En
vue d'avoir une vue d'ensemble des protéines en terme de corrélation, nous avons réalisé leurs clas-
sification hiérarchique ascendante (voir Figure 5.4) en prenant – log(p-valeur) comme dissimila-
rité. Sur la Figure 5.4, se distinguent clairement deux branches : celle de l'IgG1 et celle de l'IgG3.
En d'autres mots, les protéines IgG1_A1, IgG1_A2, IgG1_B et IgG1_C sont sur une branche dis-
tincte de celle qui porte les protéines IgG3_A1, IgG3_A2, IgG3_B et IgG3_C. Par ailleurs, IgG1
et IgG3, A1 et A2, et B et C vont ensemble. Ces résultats sont biologiquement cohérents et sont
compatibles avec la littérature, puisque A1 et A2 sont deux domaines de la même protéine, et B
et C sont des protéines différents. Sur l'arbre (Figure 5.4), il apparaît que les protéines IgG3_A1
et IgG1_B qui ne sont pas significativement corrélées d'après le test de corrélation bi-varié sont
éloignées l'une de l'autre. Statistiquement, le modèle qui sera utilisé pour l'analyse jointe des huit
protéines ne sera probablement pas celui qui contient vingt-sept corrélations significatives, évitant
un possible sur-ajustement du modèle. En s'appuyant sur ces résultats, il semble digne d'intérêt de
recourir à une procédure de régularisation dans l'ajustement du modèle conjoint à huit variables
dépendantes.

5.3.2 Deuxième illustration : Estimation par optimisation de la déviance profilée du modèle

Dans cette section, nous utilisons la procédure d'estimation introduite au Chapitre 3 pour illus-
trer l'estimation des paramètres d'un modèle bi-varié ajusté aux données et ayant pour variables
dépendantes conc.IgG1_A1 et conc.IgG3_A2 comme suit.

$$
\begin{align}
\text{conc.IgG1\_A1} &= (\mathbf{1}, \text{ap}, \text{conc\_CO.IgG1\_A1}, \text{conc\_M3.IgG1\_A1}, \text{hb}, \text{inf\_trim}, \nonumber \\
&\phantom{=} \text{pred\_trim}, \text{nutri\_trim})\beta_1 + (\mathbf{1}, \text{inf\_trim})\gamma_1 + \varepsilon_1, \nonumber \\
\text{conc.IgG3\_A2} &= (\mathbf{1}, \text{ap}, \text{conc\_CO.IgG3\_A2}, \text{conc\_M3.IgG3\_A2}, \text{hb}, \text{inf\_trim}, \nonumber \\
&\phantom{=} \text{pred\_trim}, \text{nutri\_trim})\beta_2 + (\mathbf{1}, \text{inf\_trim})\gamma_2 + \varepsilon_2, \tag{5.2}
\end{align}
$$

avec

$$
\gamma = (\gamma_1^\top, \gamma_2^\top)^\top \sim \mathcal{N} \left( \mathbf{0}, \bar{\Gamma} \right), \quad \varepsilon = (\varepsilon_1^\top, \varepsilon_2^\top)^\top \sim \mathcal{N} \left( \mathbf{0}, \begin{pmatrix} \sigma_1^2 \mathbf{I} & 0 \\ 0 & \sigma_2^2 \mathbf{I} \end{pmatrix} \right). \quad (5.3)
$$

Ici, notre stratégie a consisté à 1) estimer les paramètres du modèle en démarrant la procédure à partir de vingt-cinq points de départ aléatoirement choisis, et de 2) considérer l'estimation ayant la plus grande valeur de vraisemblance comme fournissant la "vraie" valeur du vecteur de paramètres recherché. Puis nous calculons l'erreur quadratique moyenne (MSE) des vingt-quatre estimations restantes. Ceci permettrait d'évaluer la sensibilité de la procédure d'estimation au point de départ de l'algorithme, dans un cadre d'analyse de données réelles. Les résultats obtenus sont contenus dans le tableau 5.2. A partir de ces résultats, nous constatons effectivement que la sensibilité de notre procédure d'estimation au point de départ est vraiment minime, contrairement à beaucoup d'autres procédures dans la littérature. Par ailleurs, l'estimateur de la matrice de covariance des
<table><caption>TABLE 5.2 : Analyse empirique de données de paludisme.</caption><thead><tr><th rowspan="2">Covariables</th><th colspan="4">Variables réponse</th></tr><tr><th colspan="2">conc.IgG1_A1</th><th colspan="2">conc.IgG3_A2</th></tr><tr><th></th><th>β̂<sub>1</sub></th><th>MSE</th><th>β̂<sub>2</sub></th><th>MSE</th></tr></thead><tbody><tr><td>Intercept</td><td>0.609</td><td>9.05 × 10<sup>-5</sup></td><td>-1.626</td><td>3.18 × 10<sup>-5</sup></td></tr><tr><td>ap</td><td>-0.093</td><td>1.06 × 10<sup>-5</sup></td><td>-0.337</td><td>1.04 × 10<sup>-6</sup></td></tr><tr><td>conc_CO.IgG1_A1</td><td>0.160</td><td>1.68 × 10<sup>-6</sup></td><td>—</td><td>—</td></tr><tr><td>conc_M3.IgG1_A1</td><td>0.148</td><td>9.85 × 10<sup>-6</sup></td><td>—</td><td>—</td></tr><tr><td>conc_CO.IgG3_A2</td><td>—</td><td>—</td><td>0.047</td><td>6.44 × 10<sup>-7</sup></td></tr><tr><td>conc_M3.IgG3_A2</td><td>—</td><td>—</td><td>0.155</td><td>2.22 × 10<sup>-7</sup></td></tr><tr><td>hb</td><td>-0.162</td><td>3.22 × 10<sup>-7</sup></td><td>-0.345</td><td>1.35 × 10<sup>-7</sup></td></tr><tr><td>inf_trim</td><td>0.369</td><td>1.89 × 10<sup>-6</sup></td><td>0.696</td><td>5.09 × 10<sup>-7</sup></td></tr><tr><td>pred_trim</td><td>-0.003</td><td>5.25 × 10<sup>-8</sup></td><td>0.017</td><td>1.49 × 10<sup>-8</sup></td></tr><tr><td>nutri_trim</td><td>0.024</td><td>5.81 × 10<sup>-6</sup></td><td>0.115</td><td>3.71 × 10<sup>-5</sup></td></tr></tbody><tfoot><tr><td></td><td>σ̂<sub>1</sub> = 1.395</td><td>4.96 × 10<sup>-6</sup></td><td>σ̂<sub>2</sub> = 1.626</td><td>2.42 × 10<sup>-5</sup></td></tr></tfoot></table>

effets aléatoires est donné par

$$ \hat{\Gamma} = \begin{pmatrix} 0.58 & -0.13 & 0.74 & -0.36 \\ -0.13 & 0.23 & -0.39 & 0.37 \\ 0.74 & -0.39 & 0.94 & -0.24 \\ -0.36 & 0.37 & -0.24 & 0.34 \end{pmatrix}, \quad (5.4) $$

avec une MSE de 0.0095.
# CONCLUSION

## 6.1 Bilan

La question d'analyse jointe et longitudinale de plusieurs variables statistiques est d'actualité et suscite un intérêt croissant. Dans cette thèse, nous avons pu nous y pencher dans le cadre précis d'analyse de données longitudinales multidimensionnelles relatives à l'étude d'acquisition immunitaire contre le paludisme chez l'enfant en Afrique au sud du Sahara, plus précisément au Bénin. Notre démarche a été d'explorer dans un premier temps la littérature sur le sujet, en matière de démarches méthodologiques et algorithmiques. Ce qui nous a permis d'avoir une vue d'ensemble sur les différentes méthodes existantes afin de nous positionner en terme de contributions scientifiques.

Le chapitre 1 a donc fait l'objet d'un exposé du développement historique des méthodes d'analyse de données longitudinales unidimensionnelles où l'on analyse les caractéristiques de la distribution d'une seule variable longitudinale, conditionnellement à un ou plusieurs facteurs d'intérêt. Ce chapitre nous a aussi permis de proposer une définition plus générale du modèle linéaire à effets mixtes que celle rencontrée dans la littérature.

Le chapitre 2 expose nos travaux publiés dans PLoS One et qui introduisent des expressions plus générales (que celles rencontrées dans la littérature) d'estimateurs du maximum de vraisemblance des paramètres du modèle linéaire multidimensionnel à effets mixtes, utilisant l'algorithme EM. En effet, ces estimateurs intègrent naturellement le cadre d'analyse de données multidimensionnelles multi-niveaux dont les données longitudinales multidimensionnelles peuvent être vues comme étant un cas particulier. Nous avons également introduit un test de corrélation bi-varié permettant de tester la significativité globale des corrélations entre les effets aléatoires de deux différentes dimensions du modèle. Ce test a pour but d'aider à savoir si deux variables dépendantes peuvent être analysées séparément ou non, puisque la liaison entre elles (variables dépendantes) dépend uniquement des effets aléatoires, hypothèse faite à dessein. Ceci permettra de construire un modèle conjoint plus parcimonieux en terme de composantes de variance des effets aléatoires à travers une procédure de sélection pas-à-pas ascendante. A ce niveau, bien que la distribution asymptotique de la statistique de test soit empiriquement un $\chi^2$ sur nos simulations, la littérature pointe plutôt du doigt un mélange pondéré de $\chi^2$. Ce qui nous amène à envisager de considérer un approfondissement de la question dans nos prochains travaux à travers une démarche méthodologique purement théorique, ne serait-ce que dans le cadre spécifique de ce test de corrélation bi-varié.

En ce qui concerne le chapitre 3, nous y apportons une contribution méthodologique quant à l'estimation des paramètres d'un modèle linéaire multidimensionnel à effets mixtes où les résidus dimensionnels sont homoscédastiques avec la possibilité d'avoir des variances différentes dans
différentes dimensions du modèle. La procédure proposée ici est la généralisation de celle utilisée
par les auteurs du package lme4 du logiciel R. Sa sensibilité au point de départ est quasi inexistante
à la fois sur les données simulées que sur les données d'immunité palustre. Ce qui représente l'une
des forces de la méthode qui, d'ailleurs, utilise les optimiseurs existants sous le logiciel R (où nous
avons effectué tous nos travaux) pour minimiser ou maximiser le critère choisi (ML our REML)
sans efforts supplémentaires. Les estimations réalisées sur des données simulées ont montré très
clairement, une décroissance nette vers zéro de l'erreur quadratique moyenne lorsque le nombre
d'observations augmente. Cette approche semble donc fournir des estimateurs consistants, même
pour les composantes de variance. Cet aspect (la consistance des estimateurs) mériterait d'être
théoriquement approfondi dans nos travaux à venir. Enfin, le chapitre 4 quant à lui nous a permis
de faire de la sélection d'effets fixes dans la version unidimensionnelle du modèle par pénalisation
de la vraisemblance profilée des paramètres à l'exception de la variance résiduelle. La pénalité
introduite ici est de type adaptive ridge itérative permettant d'approximer les performances d'une
pénalité de type $L_0$ de la vraisemblance. Cette procédure de sélection mériterait d'être étendue à
la version multidimensionnelle du modèle où on pourra sélectionner à la fois des effets fixes et des
effets aléatoires en grande dimension (le nombre de covariables dépassant le nombre d'observa-
tions). L'approche utilisée dans les chapitres 3 et 4 mériterait sans doute d'intégrer ultérieurement
une structure plus complexe de la matrice de covariance des résidus en présence d'une corrélation
non nulle avec les effets aléatoires tout en conservant ses performances actuelles. Ce qui fera l'ob-
jet de nos réflexions dans nos travaux à venir.

## 6.2 Perspectives

La première perspective de notre travail concerne la loi asymptotique et à distance finie de la statis-
tique du rapport de vraisemblance utilisée dans le chapitre 2. Si par exemple, $\gamma_1$ et $\gamma_2$ sont respecti-
vement les vecteurs des effets aléatoires dans les dimensions 1 et 2 du modèle, le test de corrélation
bi-varié du chapitre 2 a consisté à tester $H_0 : \text{Cor}(\gamma_1, \gamma_2) = 0$ contre $H_1 : \text{Cor}(\gamma_1, \gamma_2) \neq 0$, où
nous avons procédé à de nombreuses simulations afin d'approcher numériquement la distribution
asymptotique de la statistique de test du rapport de vraisemblance sous $H_0$. Nos simulations in-
diquent que la distribution asymptotique est celle du $\chi^2(df)$ (voir Figure 6.1). Cependant, dans la
littérature, plusieurs travaux théoriques prouvent que la statistique de test du rapport de vraisem-
blance dans des conditions non standards (voir par exemple [31, 74, 171, 196]) où le paramètre
testé est sujet à contraintes, suit asymptotiquement un mélange pondéré de $\chi^2$ sous $H_0$. Or le test
de corrélation bi-varié qui nous intéresse ici porte sur les corrélations entre les effets aléatoires de
deux différentes dimensions ; lesquelles corrélations sont bien sûr sujettes à la contrainte d'ap-
partenir à l'intervalle $[-1, 1]$. Dans la suite de nos travaux, nous nous pencherons donc plus
théoriquement sur la distribution asymptotique de la statistique de test dans le cas précis de ce
test de corrélation bi-varié.

Lorsqu'on travaille avec des jeux de données de taille très modeste, on peut recourir à une
correction de type Bartlett [12] de la statistique de test *S* du rapport de vraisemblance afin de
s'assurer, tout au moins, que la statistique utilisée (c'est à dire la statistique de test corrigée *S*<sub>B</sub>) a
df comme espérance mathématique, sous l'hypothèse *H*<sub>0</sub>. Dans le chapitre 2, nous avons pratiqué
<figure><img src="image_17.png" /><figcaption>FIGURE 6.1 : Analyse empirique de la distribution asymptotique de la statistique du rapport de vraisemblance S sous $H_0$, utilisant 3000 replications de jeux de données longitudinales de tailles $N = 15000$ (observations) et $n = 500$ (subjects), simulées sous $H_0$. Les minima et les maxima de $1000 \times 3000$ réalisations simulées du $\chi^2(4)$ sont utilisés pour construire l'enveloppe rouge. La courbe bleue représente les statistiques de rapport de vraisemblance des tests de corrélation bi-varié réalisés sur les données simulées.</figcaption></figure>

une correction empirique de Bartlett ($\hat{S}_B = df \times S/\hat{E}[S|H_0]$) pour les jeux de données de petite taille. Dans la suite de nos travaux, nous entreprendrons de calculer explicitement l'espérance conditionnelle $\hat{E}[S|H_0]$, afin d'avoir une expression théorique générale pour la statistique corrigée dans le cas du test de corrélation bi-varié.

L'estimation des paramètres (du modèle) à la lme4 fait l'hypothèse de résidus homosédastiques. Cette hypothèse, quoique restrictive, nous permet d'obtenir de très bonnes performances en ce qui concerne la qualité des estimations. Dans la suite de nos travaux, nous pensons qu'il sera très utile de chercher à lever cette restriction en fournissant une expression générale de l'estimateur de la matrice de variance-covariance $\Sigma_k$ des résidus de la dimension $k$ du modèle.

La procédure de sélection d'effets fixes proposée au chapitre 4 fera l'objet d'une extension dans la suite de nos travaux. Cette extension prendra en compte la sélection d'effets aléatoires aussi bien que la sélection d'effets fixes dans la version multidimensionnelle du modèle.

D'autres travaux à venir incluent l'introduction d'une nouvelle méthode de sélection de SNP dans le traitement des données génome wide, où la procédure *adaptive ridge itérative* proposée au chapitre 4 sera utilisée. Par ailleurs, nous comptons aussi nous intéressés à l'analyse jointe de données longitudinales et de survie avec prise en compte ou non d'événements récurrents.
Bibliographie

[1] E. H. Adjakossa, I. Sadissou, M. N. Hounkonnou, and G. Nuel. Multivariate longitudinal analysis with bivariate correlation test. *PloS one*, 11(8):e0159649, 2016.

[2] G. B. Airy. *On the algebraical and numerical theory of errors of observations and the combination of observations*. Macmillan&Company, 1861.

[3] H. Akaike. Information theory and an extension of the maximum likelihood principle. *Second International Symposium on Information Theory*, (Edited by B. N. Petrov and F. Csaki):267-281, 1973.

[4] H. Akaike. A new look at the statistical model identification. *IEEE transactions on automatic control*, 19(6):716-723, 1974.

[5] P. L. Alonso, J. Sacarlal, J. J. Aponte, A. Leach, E. Macete, J. Milman, I. Mandomando, B. Spiessens, C. Guinovart, M. Espasa, et al. Efficacy of the rts, s/as02a vaccine against plasmodium falciparum infection and disease in young african children : randomised controlled trial. *The Lancet*, 364(9443):1411-1420, 2004.

[6] X. An, Q. Yang, and P. M. Bentler. A latent factor linear mixed model for high-dimensional longitudinal data analysis. *Statistics in medicine*, 32(24):4229-4239, 2013.

[7] R. F. Anders, C. G. Adda, M. Foley, and R. S. Norton. Recombinant protein vaccines against the asexual blood-stages of plasmodium falciparum. *Human vaccines*, 6(1):39-53, 2010.

[8] R. H. Baayen, D. J. Davidson, and D. M. Bates. Mixed-effects modeling with crossed random effects for subjects and items. *Journal of memory and language*, 59(4):390-412, 2008.

[9] J. Ballard, J. Khoury, K. Wedig, L. Wang, B. Eilers-Walsman, and R. Lipp. New ballard score, expanded to include extremely premature infants. *The Journal of pediatrics*, 119(3): 417-423, 1991.

[10] S. Bandyopadhyay, B. Ganguli, and A. Chatterjee. A review of multivariate longitudinal data analysis. *Statistical methods in medical research*, 20(4):299-330, 2011.

[11] A. Barron, L. Birgé, and P. Massart. Risk bounds for model selection via penalization. *Probability theory and related fields*, 113(3):301-413, 1999.
[12] M. S. Bartlett. Properties of sufficiency and statistical tests. *Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences*, pages 268–282, 1937.

[13] D. Bates. Computational methods for mixed models. *LME4: Mixed-Effects Modeling with R*, pages 99–118, 2014.

[14] D. Bates, D. Sarkar, M. D. Bates, and L. Matrix. The lme4 package. *R package version*, 2 (1):74, 2007.

[15] D. Bates, M. Maechler, B. Bolker, and S. Walker. *lme4: Linear mixed-effects models using Eigen and S4*, 2013. URL http://CRAN.R-project.org/package=lme4. R package version 1.0-5.

[16] D. Bates, M. Maechler, B. Bolker, and S. Walker. lme4: Linear mixed-effects models using eigen and s4. r package version 1.1-7. *This is computer program (R package). The URL of the package is: http://CRAN.R-project.org/package=lme4*, 2014.

[17] D. M. Bates. lme4: Mixed-effects modeling with r. URL http://lme4.r-forge.r-project.org/book, 2010.

[18] L. Beckett, D. Tancredi, and R. Wilson. Multivariate longitudinal models for complex change processes. *Statistics in medicine*, 23(2):231–239, 2004.

[19] P. M. Bentler and D. G. Weeks. Linear structural equations with latent variables. *Psychometrika*, 45(3):289–308, 1980.

[20] R. Biihlmann and L. Meier. Discussion of "one-step sparse estimates in nonconcave penalized likelihood models," by h. zou and r. li. *Ann. Statist*, 36:1534–1541, 2008.

[21] H. Bozdogan. Model selection and akaike's information criterion (aic) : The general theory and its analytical extensions. *Psychometrika*, 52(3):345–370, 1987.

[22] H. Brandsma and J. Knuver. Effects of school and classroom characteristics on pupil progress in language and arithmetic. *International Journal of Educational Research*, 13(7):777–788, 1989.

[23] L. Breiman. Better subset regression using the nonnegative garrote. *Technometrics*, 37(4):373–384, 1995.

[24] L. Breiman et al. Heuristics of instability and stabilization in model selection. *The annals of statistics*, 24(6):2350–2383, 1996.

[25] L. F. Bringmann, N. Vissers, M. Wichers, N. Geschwind, P. Kuppens, F. Peeters, D. Borsboom, and F. Tuerlinckx. A network approach to psychopathology : new insights into clinical longitudinal data. *PloS one*, 8(4):e60188, 2013.

[26] C. Brombin, C. Di Serio, and P. M. Rancoita. Joint modeling of hiv data in multicenter observational studies : A comparison among different approaches. *Statistical methods in medical research*, page 0962280214526192, 2014.
[27] P. C. Bull, B. S. Lowe, M. Kortok, C. S. Molyneux, C. I. Newbold, and K. Marsh. Parasite antigens on the infected red cell surface are targets for naturally acquired immunity to malaria. *Nature medicine*, 4(3):358, 1998.

[28] E. J. Candes, M. B. Wakin, and S. P. Boyd. Enhancing sparsity by reweighted l1 minimization. *Journal of Fourier analysis and applications*, 14(5-6):877-905, 2008.

[29] E. Capanna. Grassi versus ross: who solved the riddle of malaria? *International Microbiology*, 9(1):69-74, 2006.

[30] V. J. Carey and B. A. Rosner. Analysis of longitudinally observed irregularly timed multivariate outcomes: regression with focus on cross-component correlation. *Statistics in medicine*, 20(1):21-31, 2001.

[31] D. Chant. On asymptotic tests of composite hypotheses in nonstandard conditions. *Biometrika*, 61(2):291-298, 1974.

[32] R. Charnigo, R. Kryscio, M. T. Bardo, D. Lynam, and R. S. Zimmerman. Joint modeling of longitudinal data in multiple behavioral change. *Evaluation & the health professions*, 34 (2):181-200, 2011.

[33] Z. Chen and D. B. Dunson. Random effects selection in linear mixed models. *Biometrics*, 59(4):762-769, 2003.

[34] S. Cohen, I. McGregor, S. Carrington, et al. Gamma-globulin and acquired immunity to human malaria. *Nature*, 192:733-7, 1961.

[35] V. Corbel, R. N ?guessan, C. Brengues, F. Chandre, L. Djogbenou, T. Martin, M. Akogbeto, J.-M. Hougard, and M. Rowland. Multiple insecticide resistance mechanisms in anopheles gambiae and culex quinquefasciatus from benin, west africa. *Acta tropica*, 101(3):207-216, 2007.

[36] P. Corran, P. Coleman, E. Riley, and C. Drakeley. Serology: a robust indicator of malaria transmission intensity? *Trends in parasitology*, 23(12):575-582, 2007.

[37] G. Cottrell, B. Kouwaye, C. Pierrat, A. Le Port, A. Bouraïma, N. Fonton, M. N. Hounkonnou, A. Massougbodji, V. Corbel, and A. Garcia. Modeling the influence of local environmental factors on malaria transmission in benin and its implications for cohort study. *PLoS One*, 7(1):e28812, 2012.

[38] D. Courtin, M. Oesterholt, H. Huismans, K. Kusi, J. Milet, C. Badaut, O. Gaye, W. Roeffen, E. J. Remarque, R. Sauerwein, et al. The quantity and quality of african children's igg responses to merozoite surface antigens reflect protection against plasmodium falciparum malaria. *PloS one*, 4(10):e7590, 2009.

[39] R. Crouchley, D. Stott, J. Pritchard, and D. Grose. Multivariate generalised linear mixed models via sabrer (sabre in r). 2010.
[40] M. Crowder. On the use of a working correlation matrix in using generalized linear models for repeated measures. *Biometrika*, 82(2):407-410, 1995.

[41] M. Danis and J. Mouchet. *Paludisme*. 1991.

[42] M. Davidian and D. M. Giltinan. Nonlinear models for repeated measurement data : an overview and update. *Journal of Agricultural, Biological, and Environmental Statistics*, 8 (4):387-419, 2003.

[43] C. Dechavanne. *Construction de la réponse anticorps spécifique du paludisme chez le jeune enfant : étude combinée de l'hôte, du parasite et de leur environnement*. PhD thesis, Paris 5, 2012.

[44] A. P. Dempster, N. M. Laird, and D. B. Rubin. Maximum likelihood from incomplete data via the em algorithm. *Journal of the royal statistical society. Series B (methodological)*, pages 1-38, 1977.

[45] P. Diggle. *Analysis of longitudinal data*. Oxford University Press, 2002.

[46] P. Diggle, P. Heagerty, K.-Y. Liang, and S. Zeger. *Analysis of longitudinal data*. Number 25. Oxford University Press, 2013.

[47] A. Djènontin, S. Bio-Bangana, N. Moiroux, M.-C. Henry, O. Bousari, J. Chabi, R. Ossè, S. Koudénoukpo, V. Corbel, M. Akogbéto, et al. Culicidae diversity, malaria transmission and insecticide resistance alleles in malaria vectors in ouidah-kpomasse-tori district from benin (west africa) : A pre-intervention study. *Parasit Vectors*, 3:83, 2010.

[48] R. Djouaka, H. Irving, Z. Tukur, and C. S. Wondji. Exploring mechanisms of multiple insecticide resistance in a population of the malaria vector anopheles funestus in benin. *PLos one*, 6(11):e27760, 2011.

[49] S. C. Duncan and T. E. Duncan. A multivariate latent growth curve analysis of adolescent substance use. *Structural Equation Modeling : A Multidisciplinary Journal*, 3(4):323-347, 1996.

[50] L. J. Edwards, K. E. Muller, R. D. Wolfinger, B. F. Qaqish, and O. Schabenberger. An r2 statistic for fixed effects in the linear mixed model. *Statistics in medicine*, 27(29):6137-6157, 2008.

[51] J. Fan. Comments on wavelets in statistics : A review by a. antoniadis. *Journal of the Italian Statistical Society*, 6(2):131-138, 1997.

[52] J. Fan and R. Li. Variable selection via nonconcave penalized likelihood and its oracle properties. *Journal of the American statistical Association*, 96(456):1348-1360, 2001.

[53] J. Fan and J. Lv. A selective overview of variable selection in high dimensional feature space. *Statistica Sinica*, 20(1):101, 2010.
[54] Y. Fan and R. Li. Variable selection in linear mixed effects models. *Annals of statistics*, 40 (4):2043, 2012.

[55] J. J. Faraway. *Extending the linear model with R : generalized linear, mixed effects and nonparametric regression models*. CRC press, 2005.

[56] S. Fieuws and G. Verbeke. Joint modelling of multivariate longitudinal profiles : pitfalls of the random-effects approach. *Statistics in Medicine*, 23(20):3093–3104, 2004.

[57] S. Fieuws and G. Verbeke. Pairwise fitting of mixed models for the joint modeling of multivariate longitudinal profiles. *Biometrics*, 62(2):424–431, 2006.

[58] S. Fieuws and G. Verbeke. Joint models for high-dimensional longitudinal data. *Longitudinal data analysis*, pages 367–391, 2009.

[59] S. Fieuws, G. Verbeke, B. Maes, and Y. Vanrenterghem. Predicting renal graft failure using multivariate longitudinal profiles. *Biostatistics*, 9(3):419–431, 2008.

[60] R. A. Fisher. The correlation between relatives on the supposition of mendelian inheritance. *Transactions of the royal society of Edinburgh*, 52(02):399–433, 1919.

[61] G. Fitzmaurice, M. Davidian, G. Verbeke, and G. Molenberghs. *Longitudinal data analysis*. CRC Press, 2008.

[62] G. M. Fitzmaurice, N. M. Laird, and J. H. Ware. *Applied longitudinal analysis*, volume 998. John Wiley & Sons, 2012.

[63] L. E. Frank and J. H. Friedman. A statistical view of some chemometrics regression tools. *Technometrics*, 35(2):109–135, 1993.

[64] E. W. Frees. *Longitudinal and panel data : analysis and applications in the social sciences*. Cambridge University Press, 2004.

[65] F. Frommlet and G. Nuel. An adaptive ridge procedure for 10 regularization. *PloS one*, 11 (2):e0148620, 2016.

[66] W. J. Fu. Penalized regressions : the bridge versus the lasso. *Journal of computational and graphical statistics*, 7(3):397–416, 1998.

[67] I. Funatogawa, T. Funatogawa, and Y. Ohashi. An autoregressive linear mixed effects model for the analysis of longitudinal data which show profiles approaching asymptotes. *Statistics in medicine*, 26(9):2113–2130, 2007.

[68] A. T. Galecki. General class of covariance structures for two or more repeated factors in longitudinal data analysis. *Communications in Statistics-Theory and Methods*, 23(11): 3105–3119, 1994.

[69] B. Gamain, S. Gratepanche, L. H. Miller, and D. I. Baruch. Molecular basis for the dicho-tomy in plasmodium falciparum adhesion to cd36 and chondroitin sulfate a. *Proceedings of the National Academy of Sciences*, 99(15):10020–10024, 2002.
[70] K. M. Gates and P. C. Molenaar. Group search algorithm recovers effective connectivity maps for individuals in homogeneous and heterogeneous samples. *Neuroimage*, 63(1):310–319, 2012.

[71] M. L. Gatton, J. M. Peters, E. V. Fowler, and Q. Cheng. Switching rates of plasmodium falciparum var genes: faster than we thought? *Trends in parasitology*, 19(5):202–208, 2003.

[72] A. Gelman and J. Hill. *Data analysis using regression and multilevel/hierarchical models*. Cambridge University Press, 2006.

[73] H. Geys, G. Molenberghs, and L. M. Ryan. Pseudolikelihood modeling of multivariate outcomes in developmental toxicology. *Journal of the American Statistical Association*, 94 (447):734–745, 1999.

[74] V. Giampaoli and J. M. Singer. Likelihood ratio tests for variance components in linear mixed models. *Journal of Statistical Planning and Inference*, 139(4):1435–1448, 2009.

[75] A. S. Goldberger. Best linear unbiased prediction in the generalized linear regression model. *Journal of the American Statistical Association*, 57(298):369–375, 1962.

[76] H. Goldstein. Multilevel statistical models halsted press. *New York*, 1995.

[77] M. F. Good and D. L. Doolan. Malaria vaccine design: immunological considerations. *Immunity*, 33(4):555–566, 2010.

[78] Y. Grandvalet. Least absolute shrinkage is equivalent to quadratic penalization. In *ICANN* 98, pages 201–206. Springer, 1998.

[79] S. M. Gray and R. Brookmeyer. Estimating a treatment effect from multidimensional longitudinal data. *Biometrics*, pages 976–988, 1998.

[80] S. M. Gray and R. Brookmeyer. Multidimensional longitudinal data: estimating a treatment effect from continuous, discrete, or time-to-event response variables. *Journal of the American Statistical Association*, 95(450):396–406, 2000.

[81] F. Gumedze and T. Dunne. Parameter estimation and inference in the linear mixed model. *Linear Algebra and its Applications*, 435(8):1920–1944, 2011.

[82] U. Halekoh, S. Højsgaard, and J. Yan. The r package geepack for generalized estimating equations. *Journal of Statistical Software*, 15(2):1–11, 2006.

[83] J. D. Hamilton. State-space models. *Handbook of econometrics*, 4:3039–3080, 1994.

[84] G. R. Hancock, W.-L. Kuo, and F. R. Lawrence. An illustration of second-order latent growth models. *Structural Equation Modeling*, 8(3):470–489, 2001.

[85] E. J. Hannan and B. G. Quinn. The determination of the order of an autoregression. *Journal of the Royal Statistical Society. Series B (Methodological)*, pages 190–195, 1979.
[86] H. O. Hartley and J. N. Rao. Maximum-likelihood estimation for the mixed analysis of variance model. *Biometrika*, 54(1-2):93-108, 1967.

[87] D. Harville. Extension of the gauss-markov theorem to include the estimation of random effects. *The Annals of Statistics*, pages 384–395, 1976.

[88] D. A. Harville. Maximum likelihood approaches to variance component estimation and to related problems. *Journal of the American Statistical Association*, 72(358):320–338, 1977.

[89] S. I. Hay, C. A. Guerra, A. J. Tatem, A. M. Noor, and R. W. Snow. The global distribution and population at risk of malaria: past, present, and future. *The Lancet Infectious Diseases*, 4(6):327–336, 2004.

[90] D. Hedeker and R. D. Gibbons. *Longitudinal data analysis*, volume 451. John Wiley & Sons, 2006.

[91] C. R. Henderson. Estimation of genetic parameters. In *Biometrics*, volume 6, pages 186–187. INTERNATIONAL BIOMETRIC SOC 1441 I ST, NW, SUITE 700, WASHINGTON, DC 20005-2210, 1950.

[92] C. R. Henderson. Estimation of variance and covariance components. *Biometrics*, 9(2):226–252, 1953.

[93] C. Horváth and J. E. Wieringa. Pooling data for the analysis of dynamic marketing systems. *Statistica Neerlandica*, 62(2):208–229, 2008.

[94] J. G. Ibrahim, H. Zhu, and N. Tang. Model selection criteria for missing-data problems using the em algorithm. *Journal of the American Statistical Association*, 2008.

[95] J. G. Ibrahim, H. Zhu, R. I. Garcia, and R. Guo. Fixed and random effects selection in mixed effects models. *Biometrics*, 67(2):495–503, 2011.

[96] R. I. Jennrich and M. D. Schluchter. Unbalanced repeated-measures models with structured covariance matrices. *Biometrics*, pages 805–820, 1986.

[97] K. L. Jensen, H. Spiild, and J. Toftum. Implementation of multivariate linear mixed-effects models in the analysis of indoor climate performance experiments. *International Journal of Biometeorology*, 56(1):129–136, 2012.

[98] J. Jiang, J. S. Rao, Z. Gu, T. Nguyen, et al. Fence methods for mixed model selection. *The Annals of Statistics*, 36(4):1669–1692, 2008.

[99] R. A. Johnson, D. W. Wichern, and P. Education. *Applied multivariate statistical analysis*, volume 4. Prentice hall Englewood Cliffs, NJ, 2007.

[100] R. Kobbe, R. Neuhoff, F. Marks, S. Adjei, I. Langefeld, C. Von Reden, O. Adjei, C. G. Meyer, and J. May. Seasonal variation and high multiplicity of first plasmodium falciparum infections in children from a holoendemic area in ghana, west africa. *Tropical Medicine & International Health*, 11(5):613–619, 2006.
[101] M. Kramer. R 2 statistics for mixed models. 2005.

[102] E. Kuhn and M. Lavielle. Maximum likelihood estimation in nonlinear mixed effects models. *Computational Statistics & Data Analysis*, 49(4):1020-1038, 2005.

[103] N. Laird, N. Lange, and D. Stram. Maximum likelihood computations with repeated measures: application of the em algorithm. *Journal of the American Statistical Association*, 82 (397):97-105, 1987.

[104] N. M. Laird. Computation of variance components using the em algorithm. *Journal of Statistical Computation and Simulation*, 14(3-4):295-303, 1982.

[105] N. M. Laird and J. H. Ware. Random-effects models for longitudinal data. *Biometrics*, pages 963-974, 1982.

[106] P. Lambert and F. Vandenhende. A copula-based model for multivariate non-normal longitudinal data: analysis of a dose titration safety study on a new antidepressant. *Statistics in medicine*, 21(21):3197-3217, 2002.

[107] P. Lazarsfeld and M. Fiske. The "panel" as a new tool for measuring opinion. *The Public Opinion Quarterly*, 2(4):596-612, 1938.

[108] A. Le Port, G. Cottrell, Y. Martin-Prevel, F. Migot-Nabias, M. Cot, and A. Garcia. First malaria infections in a cohort of infants in benin: biological, environmental and genetic determinants. description of the study site, population methods and preliminary results. *BMJ open*, 2(2):e000342, 2012.

[109] Y. Lee and J. Nelder. Generalized linear models for the analysis of quality-improvement experiments. *Canadian Journal of Statistics*, 26(1):95-105, 1998.

[110] K.-Y. Liang and S. L. Zeger. Longitudinal data analysis using generalized linear models. *Biometrika*, 73(1):13-22, 1986.

[111] K.-Y. Liang, S. L. Zeger, and B. Qaqish. Multivariate regression analyses for categorical data. *Journal of the Royal Statistical Society. Series B (Methodological)*, pages 3-40, 1992.

[112] X. Lin. Variance component testing in generalised linear models with random effects. *Biometrika*, 84(2):309-326, 1997.

[113] D. V. Lindley and A. F. Smith. Bayes estimates for the linear model. *Journal of the Royal Statistical Society. Series B (Methodological)*, pages 1-41, 1972.

[114] M. J. Lindstrom and D. M. Bates. Newton raphson and em algorithms for linear mixed-effects models for repeated-measures data. *Journal of the American Statistical Association*, 83(404):1014-1022, 1988.

[115] M. J. Lindstrom and D. M. Bates. Nonlinear mixed effects models for repeated measures data. *Biometrics*, pages 673-687, 1990.
[116] R. Littell, G. Milliken, W. Stroup, R. Wolfinger, and O. Schabenberger. Random coefficient models. *SAS system for mixed models*. Cary, NC : SAS Institute Inc, pages 253–66, 1996.

[117] R. C. Littell, G. A. Milliken, W. W. Stroup, R. D. Wolfinger, and O. Schabenberger. Sas system for mixed models cary. *Nc : sas institute*, 1996.

[118] W. Liu, Y. Li, G. H. Learn, R. S. Rudicell, J. D. Robertson, B. F. Keele, J.-B. N. Numba, C. M. Sanz, D. B. Morgan, S. Locatelli, et al. Origin of the human malaria parasite plasmodium falciparum in gorillas. *Nature*, 467(7314):420–425, 2010.

[119] T. Lodewyckx, F. Tuerlinckx, P. Kuppens, N. B. Allen, and L. Sheeber. A hierarchical state space approach to affective dynamics. *Journal of mathematical psychology*, 55(1):68–83, 2011.

[120] R. C. MacCallum, C. Kim, W. B. Malarkey, and J. K. Kiecolt-Glaser. Studying multivariate change using multilevel models and latent curve models. *Multivariate Behavioral Research*, 32(3):215–253, 1997.

[121] T. E. MaCurdy. The use of time series processes to model the error structure of earnings in a longitudinal data analysis. *Journal of econometrics*, 18(1):83–114, 1982.

[122] C. L. Mallows. Some comments on c p. *Technometrics*, 15(4):661–675, 1973.

[123] J. J. McArdle. Dynamic but structural equation modeling of repeated measures data. In *Handbook of multivariate experimental psychology*, pages 561–614. Springer, 1988.

[124] J. S. McCarthy, J. Marjason, S. Elliott, P. Fahey, G. Bang, E. Malkin, E. Tierney, H. Aked-Hurditch, C. Adda, N. Cross, et al. A phase 1 trial of msp2-c1, a blood-stage malaria vaccine containing 2 isoforms of msp2 formulated with montanide (r) isla 720. *PLoS One*, 6(9):e24413, 2011.

[125] X.-L. Meng and D. B. Rubin. Maximum likelihood estimation via the ecm algorithm : A general framework. *Biometrika*, 80(2):267–278, 1993.

[126] J. J. Miller. Asymptotic properties and computation of maximum likelihood estimates in the mixed model of the analysis of variance. Technical report, DTIC Document, 1973.

[127] R. Miller, S. Ikram, G. Armelagos, R. Walker, W. Harer, C. Shiff, D. Baggett, M. Carrigan, and S. Maret. Diagnosis of plasmodium falciparum infections in mummies using the rapid manual parasight ?-f test. *Transactions of the Royal Society of Tropical Medicine and Hygiene*, 88(1):31–32, 1994.

[128] G. Minois. *Histoire de l'avenir : des prophètes à la prospective*. Fayard, 1996.

[129] G. Molenberghs and G. Verbeke. Models for discrete longitudinal data. 2005.

[130] J. Mouchet. *Biodiversité du paludisme dans le monde*. John Libbey Eurotext, 2004.

[131] S. Müller, J. L. Scealy, A. H. Welsh, et al. Model selection in linear mixed models. *Statistical Science*, 28(2):135–167, 2013.
[132] R. B. Nelsen. *An introduction to copulas*. Springer, 1999.

[133] M. Nerlove. An essay on the history of panel data econometrics. In *Proceedings of Ninth International Conference on Panel Data, Geneva, Switzerland*, page 13, 2000.

[134] P. N. Oak. World vedic heritage : a history of histories : presenting a unique unified field theory of history that from the beginning of time the world practised vedic culture and spoke sanskrit. 1984.

[135] L. M. O'Brien and G. M. Fitzmaurice. Analysis of longitudinal multiple-source binary data using generalized estimating equations. *Journal of the Royal Statistical Society : Series C (Applied Statistics)*, 53(1):177-193, 2004.

[136] P. E. Okello, W. Van Bortel, A. M. Byaruhanga, A. Correwyn, P. Roelants, A. Talisuna, U. d ?Alessandro, and M. Coosemans. Variation in malaria transmission intensity in seven sites throughout uganda. *The American journal of tropical medicine and hygiene*, 75(2): 219-225, 2006.

[137] F. J. Oort. Three-mode models for multivariate longitudinal data. *British Journal of Mathematical and Statistical Psychology*, 54(1):49-78, 2001.

[138] W. H. Organization et al. World malaria report 2012. 2012. *Geneva : World Health Organization Google Scholar*, 2014.

[139] H. D. Patterson and R. Thompson. Recovery of inter-block information when block sizes are unequal. *Biometrika*, 58(3):545-554, 1971.

[140] J. Pinheiro and D. Bates. *Mixed-effects models in S and S-PLUS*. Springer Science & Business Media, 2006.

[141] J. Pinheiro, D. Bates, S. DebRoy, D. Sarkar, et al. Linear and nonlinear mixed effects models. *R package version*, 3:57, 2007.

[142] J. Pinheiro, D. Bates, S. DebRoy, and D. Sarkar. R core team (2014). nlme : linear and nonlinear mixed effects models. r package version 3.1-117. URL : http://cran.r-project.org/web/packages/nlme/index.html, 2014.

[143] R. L. Prentice and L. P. Zhao. Estimating equations for parameters in means and covariances of multivariate discrete and continuous responses. *Biometrics*, pages 825-839, 1991.

[144] W. Pu and X.-F. Niu. Selecting mixed-effects models based on a generalized information criterion. *Journal of multivariate analysis*, 97(3):733-758, 2006.

[145] R Core Team. *R : A Language and Environment for Statistical Computing*. R Foundation for Statistical Computing, Vienna, Austria, 2014. URL : http://www.R-project.org/.

[146] R Core Team. *R : A Language and Environment for Statistical Computing*. R Foundation for Statistical Computing, Vienna, Austria, 2015. URL : http://www.R-project.org/.
[147] R Core Team. *R: A Language and Environment for Statistical Computing*. R Foundation for Statistical Computing, Vienna, Austria, 2015. URL http://www.R-project.org/.

[148] J. Ramsay and B. Silverman. Functional data analysis. 1997, 1997.

[149] C. R. Rao. Estimation of variance and covariance components under minimum theory. *Journal of multivariate analysis*, 1(3):257-275, 1971.

[150] R. Rao and Y. Wu. A strongly consistent procedure for model selection in a regression problem. *Biometrika*, 76(2):369-374, 1989.

[151] G. Reinsel. Estimation and prediction in a multivariate random effects generalized linear model. *Journal of the American Statistical Association*, 79(386):406-414, 1984.

[152] H. Ribaudo and S. Thompson. The analysis of repeated multivariate binary quality of life data: a hierarchical model approach. *Statistical methods in medical research*, 11(1):69-83, 2002.

[153] J. A. Rice and C. O. Wu. Nonparametric mixed effects models for unequally sampled noisy curves. *Biometrics*, 57(1):253-259, 2001.

[154] R. C. Rippe, J. J. Meulman, and P. H. Eilers. Visualization of genomic changes by segmented smoothing using an l0 penalty. *PloS one*, 7(6):e38230, 2012.

[155] G. K. Robinson. That blup is a good thing: the estimation of random effects. *Statistical science*, pages 15-32, 1991.

[156] J. Rochon. Analyzing bivariate repeated measures for discrete and continuous outcome variables. *Biometrics*, pages 740-750, 1996.

[157] C. Rogier, E. Orlandi-Pradines, T. Fusai, B. Pradines, S. Briolant, and L. Almeras. [malaria vaccines: prospects and reality]. *Medecine et maladies infectieuses*, 36(8):414-422, 2006.

[158] A. Sabchareon, T. Burnouf, D. Ouattara, P. Attanath, H. Bouharoun-Tayoun, P. Chantavanich, C. Foucault, T. Chongsuphajaisiddhi, and P. Druilhe. Parasitologic and clinical human response to immunoglobulin administration in falciparum malaria. *The American journal of tropical medicine and hygiene*, 45(3):297-308, 1991.

[159] M. Sammel, X. Lin, and L. Ryan. Multivariate linear mixed models for multiple outcomes. *Statistics in medicine*, 18(17-18):2479-2492, 1999.

[160] B. R. Saville and A. H. Herring. Testing random effects in the linear mixed model using approximate bayes factors. *Biometrics*, 65(2):369-376, 2009.

[161] J. L. Schafer and R. M. Yucel. Computational strategies for multivariate linear mixed-effects models with missing values. *Journal of computational and Graphical Statistics*, 11 (2):437-457, 2002.

[162] H. Scheffe. Alternative models for the analysis of variance. *The Annals of Mathematical Statistics*, pages 251-271, 1956.
[163] J. Schelldorfer. *Immlasso*: Linear mixed-effects models with lasso. *R package version* 0.1-2, 2011.

[164] J. Schelldorfer, P. Bühlmann, G. DE, and S. VAN. Estimation for high-dimensional linear mixed-effects models using 11-penalization. *Scandinavian Journal of Statistics*, 38(2):197–214, 2011.

[165] J. Schelldorfer, L. Meier, and P. Bühlmann. *Glmmlasso*: an algorithm for high-dimensional generalized linear mixed models using ?1-penalization. *Journal of Computational and Graphical Statistics*, 23(2):460–477, 2014.

[166] P. Schlagenhauf and C. Hatz. *Traitement de secours antipaludique*: Actualisation 1997. *Médecine et hygiène*, 55(2165):1126–1127, 1997.

[167] G. Schwarz et al. Estimating the dimension of a model. *The annals of statistics*, 6(2): 461–464, 1978.

[168] S. Searle, G. Casella, and C. McCulloch. Variance components john wiley and sons. *New York, New York, USA*, 1992.

[169] S. R. Searle. An overview of variance component estimation. *Metrika*, 42(1):215–230, 1995.

[170] S. R. S. R. Searle. Linear models. Technical report, 1971.

[171] S. G. Self and K.-Y. Liang. Asymptotic properties of maximum likelihood estimators and likelihood ratio tests under nonstandard conditions. *Journal of the American Statistical Association*, 82(398):605–610, 1987.

[172] A. Shah, N. Laird, and D. Schoenfeld. A random-effects model for multiple characteristics with possibly missing data. *Journal of the American Statistical Association*, 92(438):775–779, 1997.

[173] J. Shao. An asymptotic theory for linear model selection. *Statistica Sinica*, pages 221–242, 1997.

[174] N. W. Shock, r. C. Greulich, P. T. Costa, R. Andres, E. G. Lakatta, D. Arenberg, and J. D. Tobin. Normal human aging: The baltimore longitudinal study of aging. 1984.

[175] M. Sklar. *Fonctions de répartition à n dimensions et leurs marges*. Université Paris 8, 1959.

[176] T. A. Smith, R. Leuenberger, and C. Lengeler. Child mortality and malaria transmission intensity in afrika. *Trends in parasitology*, 17(3):145–149, 2001.

[177] T. A. Snijders. *Multilevel analysis*. Springer, 2011.

[178] T. A. Snijders and R. J. Bosker. Modeled variance in two-level models. *Sociological methods & research*, 22(3):342–363, 1994.
[179] G. Snounou, A. C. Grüner, C. D. Müller-Graf, D. Mazier, and L. Rénia. The plasmodium sporozoite survives rts, s vaccination. *Trends in parasitology*, 21(10):456-461, 2005.

[180] R. W. Snow, C. A. Guerra, A. M. Noor, H. Y. Myint, and S. I. Hay. The global distribution of clinical episodes of plasmodium falciparum malaria. *Nature*, 434(7030):214-217, 2005.

[181] S. Sturtz, U. Ligges, and A. Gelman. R2winbugs: A package for running winbugs from r. *Journal of Statistical Software*, 12(3):1-16, 2005. URL http://www.jstatsoft.org.

[182] S. Subramanian, D. Kim, and I. Kawachi. Covariation in the socioeconomic determinants of self rated health and happiness : a multivariate multilevel analysis of individuals and communities in the usa. *Journal of Epidemiology and Community Health*, 59(8):664-669, 2005.

[183] J. Sy, J. Taylor, and W. Cumberland. A stochastic model for the analysis of bivariate longitudinal aids data. *Biometrics*, pages 542-555, 1997.

[184] H. Theil and A. S. Goldberger. On pure and mixed statistical estimation in economics. *International Economic Review*, 2(1):65-78, 1961.

[185] R. Thiébaut, H. Jacqmin-Gadda, G. Chêne, C. Leport, and D. Commenges. Bivariate linear mixed models using sas proc mixed. *Computer methods and programs in biomedicine*, 69 (3):249-256, 2002.

[186] R. Tibshirani. Regression shrinkage and selection via the lasso. *Journal of the Royal Statistical Society. Series B (Methodological)*, pages 267-288, 1996.

[187] R. S. Tsay. *Multivariate Time Series Analysis: With R and Financial Applications*. John Wiley & Sons, 2013.

[188] W. Tschacher and F. Ramseyer. Modeling psychotherapy process by time-series panel analysis (tspa). *Psychotherapy Research*, 19(4-5):469-481, 2009.

[189] W. Tschacher, P. Zorn, and F. Ramseyer. Change mechanisms of schema-centered group psychotherapy with personality disorder patients. *PloS one*, 7(6):e39687, 2012.

[190] A. Tseloni and C. Zarafonitou. Fear of crime and victimization a multivariate multilevel analysis of competing measurements. *European Journal of Criminology*, 5(4):387-409, 2008.

[191] F. Vaida and S. Blanchard. Conditional akaike information for mixed-effects models. *Bio-metrika*, 92(2):351-370, 2005.

[192] G. Verbeke and G. Molenberghs. Linear mixed models in practice : a sas oriented approach. 1997.

[193] G. Verbeke and G. Molenberghs. *Linear mixed models for longitudinal data*. Springer Science & Business Media, 2009.
[194] G. Verbeke, S. Fieuws, G. Molenberghs, and M. Davidian. The analysis of multivariate longitudinal data : A review. *Statistical methods in medical research*, 23(1) :42–59, 2014.

[195] E. F. Vonesh. A note on the use of laplace's approximation for nonlinear mixed-effects models. *Biometrika*, 83(2) :447–452, 1996.

[196] H. Vu, S. Zhou, et al. Generalization of likelihood ratio tests under nonstandard conditions. *The Annals of Statistics*, 25(2) :897–916, 1997.

[197] W.-L. Wang and T.-H. Fan. Ecm-based maximum likelihood inference for multivariate linear mixed models with autoregressive errors. *Computational Statistics & Data Analysis*, 54(5) :1328–1341, 2010.

[198] X.-F. Wang. Joint generalized models for multidimensional outcomes : A case study of neuroscience data from multimodalities. *Biometrical Journal*, 54(2) :264–280, 2012.

[199] G. E. Weiss, B. Traore, K. Kayentao, A. Ongoiba, S. Doumbo, D. Doumtabe, Y. Kone, S. Dia, A. Guindo, A. Traore, et al. The plasmodium falciparum-specific human memory b cell compartment expands gradually with repeated malaria infections. *PLoS Pathog*, 6(5) : e1000912, 2010.

[200] S. S. Wilks. The large-sample distribution of the likelihood ratio for testing composite hypotheses. *The Annals of Mathematical Statistics*, 9(1) :60–62, 1938.

[201] H. Wu and J.-T. Zhang. *Nonparametric regression methods for longitudinal data analysis : mixed-effects modeling approaches*, volume 515. John Wiley & Sons, 2006.

[202] R. Xu. Measuring explained variation in linear mixed effects models. *Statistics in medicine*, 22(22) :3527–3541, 2003.

[203] R. Yucel. *mlmmm : ML estimation under multivariate linear mixed models with missing values*, 2010. URL https://CRAN.R-project.org/package=mlmmm. R package version 0.3-1.2.

[204] S. Zeger, K. Liang, and P. Albert. Models for longitudinal data : a generalized estimating equation approach. *Biometrics*, pages 1049–1060, 1988.

[205] S. L. Zeger and K.-Y. Liang. Longitudinal data analysis for discrete and continuous outcomes. *Biometrics*, pages 121–130, 1986.

[206] A. Zellner. An efficient method of estimating seemingly unrelated regressions and tests for aggregation bias. *Journal of the American statistical Association*, 57(298) :348–368, 1962.

[207] C. H. Zhang. Penalized linear unbiased selection. *Department of Statistics and Bioinformatics, Rutgers University*, pages 2007–003, 2007.

[208] M. Zhang, A. A. Tsiatis, M. Davidian, K. S. Pieper, and K. W. Mahaffey. Inference on treatment effects from a randomized clinical trial in the presence of premature treatment discontinuation : the synergy trial. *Biostatistics*, 12(2) :258–269, 2011.
[209] H. Zou. The adaptive lasso and its oracle properties. *Journal of the American statistical association*, 101(476):1418–1429, 2006.

[210] A. Zuur, E. N. Ieno, N. Walker, A. A. Saveliev, and G. M. Smith. *Mixed effects models and extensions in ecology with R*. Springer Science & Business Media, 2009.
Résumé

Dans cette thèse, nous nous sommes focalisés sur le modèle statistique linéaire à effets mixtes. Nous nous sommes d'abord intéressés à l'estimation consistante des paramètres du modèle dans sa version multidimensionnelle, puis à de la sélection d'effets fixes en dimension un. En ce qui concerne l'estimation des paramètres du modèle linéaire à effets mixtes multidimensionnel, nous avons proposé des estimateurs du maximum de vraisemblance par utilisation de l'algorithme EM, mais avec des expressions plus générales que celles de la littérature classique, permettant d'analyser non seulement des données longitudinales multivariées mais aussi des données multidimensionnelles multi-niveaux. Ici, en s'appuyant sur ces EM-estimateurs, nous avons introduit un test de rapport de vraisemblance permettant de tester la significativité globale des corrélations entre les effets aléatoires de deux dimensions du modèle. Ce qui permettrait de construire un modèle multidimensionnel plus parcimonieux en terme de paramètres de variance des effets aléatoires, par une procédure de selection pas-à-pas ascendante. Cette démarche a été suscitée par le fait que la dimension du vecteur de tous les effets aléatoires du modèle peut très rapidement croître avec le nombre de variables à analyser, entrainant facilement des problèmes numériques dans l'optimisation du critère choisi (ML ou REML). Nous avons ensuite proposé une procédure d'estimation consistante des paramètres du modèle qui passe par la résolution d'un problème de moindres carrés pénalisés pour fournir une expression explicite de la déviance à minimiser. La procédure de sélection d'effets fixes proposée ici est de type adaptive ridge itérative et permet d'approximer les performances de sélection d'une pénalité de type $L_0$ de la vraisemblance des paramètres du modèle. Nos résultats ont été appuyés par des études de simulation à plusieurs niveaux, mais aussi par l'analyse de plusieurs jeux de données réelles.

Abstract

This thesis focuses on the statistical linear mixed-effects model, where we have been interested in its multivariate version's parameters estimation but also in the unidimensional selection of fixed effects. Concerning the parameters estimation of the multivariate linear mixed-effects model, we have first introduced more general expressions of the EM algorithm-based estimators which fit the multivariate longitudinal data analysis framework but also the framework of the multivariate multilevel data analysis. Since the dimensionality of the total vector of random effects in the multivariate model can grow with the number of the outcome variables leading often to computational problems in the likelihood optimization, we introduced a likelihood ratio test for testing the global effect of the correlations between the random effects of two dimensions of the model. This bivariate correlation test is intended to help in constructing a more parsimonious model regarding the variance components of the random effects, using a stepwise procedure. Secondly, we have introduced another estimation procedure that yields to consistent estimates for all the model parameters. This procedure is based on the Cholesky factorization of the random effects covariance matrix and the resolution of a preliminary penalized means square problem, and leads to an explicit expression of the profiled deviance of the model. For selecting fixed effects in the one dimensional mixed-effects model, we introduce an iterative adaptive ridge procedure for approximating
$L_0$ penalty selection performances. All the results in this manuscript have been accompanied by extensive simulation studies along with real data analysis examples.
<figure><img src="image_18.png" /></figure>