Profiled deviance for the multivariate linear mixed-effects model fitting

TL;DR

This paper introduces a novel approach for estimating the covariance matrix of random effects in multivariate linear mixed-effects models, improving accuracy over traditional methods like EM, demonstrated through simulations and real-world data analysis.

Contribution

It proposes a new profiling deviance method for better covariance estimation in multivariate mixed-effects models, addressing limitations of existing estimation techniques.

Findings

The new method provides more accurate covariance estimates than EM.

Simulation studies show improved estimation performance.

Application to real data demonstrates practical utility.

Abstract

This paper focuses on the multivariate linear mixed-effects model, including all the correlations between the random effects when the marginal residual terms are assumed uncorrelated and homoscedastic with possibly different standard deviations. The random effects covariance matrix is Cholesky factorized to directly estimate the variance components of these random effects. This strategy enables a consistent estimate of the random effects covariance matrix which, generally, has a poor estimate when it is grossly (or directly) estimated, using the estimating methods such as the EM algorithm. By using simulated data sets, we compare the estimates based on the present method with the EM algorithm-based estimates. We provide an illustration by using the real-life data concerning the study of the child's immune against malaria in Benin (West Africa).

Tables

Table 1. Table 1: Mean Square Error of estimators with 95 % percent 95 \displaystyle 95\% CI estimated on 100 100 \displaystyle 100 replications for values of n ∈ { 50 , 60 , 100 , 300 } 𝑛 50 60 100 300 \displaystyle n\in\{50,60,100,300\} and N ∈ { 600 , 1000 , 3000 } 𝑁 600 1000 3000 \displaystyle N\in\{600,1000,3000\} .

Parameter	$n$	$N=600$	$N=1000$	$N=3000$
${\beta }_1$	50	2.43 (0.11 - 7.11)	1.89 (0.22 - 4.89)	1.02 (0.07 - 2.44)
	60	2.57 (0.14 - 7.87)	2.13 (0.26 - 5.55)	0.77 (0.10 - 2.09)
	100	2.27 (0.16 - 5.61)	1.55 (0.14 - 5.17)	0.71 (0.04 - 1.85)
	300	3.16 (0.14 - 10.54)	1.70 (0.10 - 4.55)	0.51 (0.04 - 1.36)
${\beta }_2$	50	5.50(0.09 - 16.35)	3.26 (0.02 - 11.64)	2.06 (0.09 - 5.98)
	60	5.06 (0.12 - 15.09)	3.22 (0.02 - 10.24)	1.78 (0.10 - 5.62)
	100	4.33 (0.03 - 13.17)	2.37 (0.02 - 6.89)	1.06 (0.02 - 3.60)
	300	4.58 (0.18 - 15.92)	2.43 (0.05 - 7.39)	0.90 (0.04 - 2.88)
${\sigma }_1$	50	0.03 (0.00 - 0.14)	0.02 (0.00 - 0.09)	0.00 (0.00 - 0.02)
	60	0.04 (0.00 - 0.11)	0.02 (0.00 - 0.08)	0.00 (0.00 - 0.02)
	100	0.03 (0.00 - 0.13)	0.01 (0.00 - 0.06)	0.00 (0.00 - 0.01)
	300	0.05 (0.00 - 0.23)	0.03 (0.00 - 0.13)	0.00 (0.00 - 0.02)
${\sigma }_2$	50	0.04 (0.00 - 0.15)	0.03 (0.00 - 0.08)	0.00 (0.00 - 0.02)
	60	0.04 (0.00 - 0.18)	0.03 (0.00 - 0.12)	0.01 (0.00 - 0.03)
	100	0.06 (0.00 - 0.18)	0.03 (0.00 - 0.11)	0.01 (0.00 - 0.03)
	300	0.13 (0.00 - 0.45)	0.04 (0.00 - 0.15)	0.01 (0.00 - 0.05)
$\overline{𝚪}$	50	0.90 (0.12 - 2.41)	0.62 (0.06 - 1.41)	0.45 (0.04 - 1.14)
	60	1.07 (0.06 - 2.64)	0.68 (0.08 - 1.98)	0.25 (0.03 - 0.66)
	100	0.90 (0.05 - 2.40)	0.47 (0.03 - 1.29)	0.21 (0.02 - 0.71)
	300	1.45 (0.19 - 4.39)	0.69 (0.06 - 1.94)	0.17 (0.02 - 0.57)

Table 2. Table 2: Mean Square Error of estimators with 95 % percent 95 \displaystyle 95\% CI estimated on 100 100 \displaystyle 100 replications for values of ( n , N ) ∈ { ( 500 , 7000 ) , ( 600 , 8000 ) , ( 800 , 10000 ) , ( 1000 , 15000 ) } 𝑛 𝑁 500 7000 600 8000 800 10000 1000 15000 \displaystyle(n,N)\in\{(500,7000),(600,8000),(800,10000),(1000,15000)\} .

Parameter	$n=500,N=7000$	$n=600,N=8000$	$n=800,N=10000$	$n=1000,N=15000$
${\beta }_1$	0.22 (0.01 - 0.62)	0.17 (0.01 - 0.48)	0.16 (0.01 - 0.47)	0.11 (0.00 - 0.32)
${\beta }_2$	0.32 (0.00 - 1.01)	0.34 (0.01 - 1.16)	0.25 (0.02 - 0.67)	0.21 (0.00 - 0.69)
${\sigma }_1$	0.00 (0.00 - 0.00)	0.00 (0.00 - 0.01)	0.00 (0.00 - 0.00)	0.00 (0.00 - 0.00)
${\sigma }_2$	0.00 (0.00 - 0.01)	0.00 (0.00 - 0.01)	0.00 (0.00 - 0.01)	0.00 (0.00 - 0.00)
$\overline{𝚪}$	0.09 (0.00 - 0.25)	0.07 (0.00 - 0.19)	0.06 (0.00 - 0.19)	0.03 (0.00 - 0.09)

Table 3. Table 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.

		Naive initialization				Advised initialization
		Cmlme		EM		Cmlme		EM
Parameter	Value	Emp. Mean	95% CI	Emp. Mean	95% CI	Emp. Mean	95% CI	Emp. Mean	95% CI
${\beta }_1$	50.67	50.79	$49.14-52.11$	13.47	$-76.70-43.37$	50.80	$49.15-52.12$	50.78	$49.09-52.01$
	-4.80	$-5.00$	$-8.39--3.66$	-4.79	$-8.08--3.42$	-5.02	$-8.39--3.66$	-4.98	$-8.38--3.65$
	14.00	14.02	$13.27-14.45$	-2.05	$-4.70--0.40$	14.02	$13.28-14.45$	14.02	$13.27-14.45$
	2.70	2.70	$2.66-2.72$	2.69	$2.66-2.72$	2.70	$2.66-2.72$	2.70	$2.66-2.72$
${\beta }_2$	13.20	13.65	$11.79-15.06$	-84.47	$-114.28--50.63$	13.65	$11.79-15.07$	13.68	$11.81-15.14$
	-2.80	-2.80	$-4.74--0.43$	-2.75	$-4.90-0.21$	-2.81	$-4.79--0.43$	-2.85	$-4.80--0.51$
	27.00	27.00	$26.87-27.10$	0.90	$-1.62-2.68$	27.00	$26.87-27.10$	27.00	$26.87-27.10$
	1.70	1.68	$1.64-1.71$	1.68	$1.64-1.71$	1.68	$1.64-1.71$	1.68	$1.64-1.71$
${\sigma }_1$	5.80	5.79	$5.62-5.92$	5.78	$5.64-5.98$	5.78	$5.64-5.92$	5.79	$5.65-5.94$
${\sigma }_2$	7.60	7.61	$7.34-7.74$	7.59	$7.33-7.73$	7.61	$7.34-7.73$	7.63	$7.36-7.73$
Nbr. of iteration	Min	56	-	63	-	48	-	14	-
	Mean	71	-	109	-	64	-	169	-
	Max	103	-	157	-	89	-	645	-

Table 4. Table 4: EM estimates compared with Cmlme estimates on the same data sets. Empirical relative error of the estimates with their 95% CI

		Naive initialization				Advised initialization
		Cmlme		EM		Cmlme		EM
Parameter	Value	R. Error	95% CI	R. Error	95% CI	R. Error	95% CI	R. Error	95% CI
${\beta }_1$	50.67	0.01	$0.00-0.03$	0.73	$0.01-2.18$	0.01	$0.00-0.03$	0.01	$0.00-0.03$
	-4.80	$0.21$	$0.02-0.32$	0.21	$0.00-0.32$	0.21	$0.02-0.32$	0.21	$0.00-0.33$
	14.00	0.02	$0.00-0.04$	1.14	$1.01-1.32$	0.02	$0.00-0.04$	0.02	$0.00-0.04$
	2.70	0.00	$0.00-0.01$	0.00	$0.00-0.01$	0.00	$0.00-0.01$	0.00	$0.00-0.01$
${\beta }_2$	13.20	0.07	$0.00-0.14$	7.39	$4.19-9.40$	0.07	$0.00-0.14$	0.07	$0.00-0.14$
	-2.80	0.43	$0.00-0.84$	0.43	$0.02-1.07$	0.43	$0.00-0.84$	0.43	$0.00-0.81$
	27.00	0.00	$0.00-0.00$	0.96	$0.88-1.02$	0.00	$0.00-0.00$	0.00	$0.00-0.00$
	1.70	0.01	$0.00-0.02$	0.01	$0.00-0.02$	0.01	$0.00-0.02$	0.01	$0.00-0.03$
${\sigma }_1$	5.80	0.01	$0.00-0.02$	0.01	$0.00-0.03$	0.01	$0.00-0.02$	0.01	$0.00-0.03$
${\sigma }_2$	7.60	0.01	$0.00-0.02$	0.01	$0.00-0.04$	0.01	$0.00-0.02$	0.01	$0.00-0.02$

Table 5. Table 5: Variables present in the empirical dataset

N∘	Variable	Description
1	id	Child ID
2	conc.Y	concentration of Y
3	conc$\_$CO.Y	Measured concentration of Y in the umbilical cord blood
4	conc$\_$M3.Y	Predicted concentration of Y in the child’s peripheral blood at 3 months
5	ap	Placental apposition
6	hb	Hemoglobin level
7	inf$\_$trim	Number of malaria infections in the previous 3 months
8	pred$\_$trim	Quarterly average number of mosquitoes child is exposed to
9	nutri$\_$trim	Quarterly average nutrition scores

Table 6. Table 6: Empirical data analysis.

	Response variables
	$\text{conc.IgG}1\_\text{A}1$		$\text{conc.IgG}3\_\text{A}2$
Covariates	Estimate	MSE	Estimate	MSE
Intercept	$0.609$	$9.05\times {10}^{-5}$	$-1.626$	$3.18\times {10}^{-5}$
ap	$-0.093$	$1.06\times {10}^{-5}$	$-0.337$	$1.04\times {10}^{-6}$
$\text{conc}\_\text{CO.IgG1}\_\text{A}1$	$0.160$	$1.68\times {10}^{-6}$	$-$	$-$
$\text{conc}\_\text{M3.IgG1}\_\text{A}1$	$0.148$	$9.85\times {10}^{-6}$	$-$	$-$
$\text{conc}\_\text{CO.IgG3}\_\text{A}2$	$-$	$-$	$0.047$	$6.44\times {10}^{-7}$
$\text{conc}\_\text{M3.IgG3}\_\text{A}2$	$-$	$-$	$0.155$	$2.22\times {10}^{-7}$
hb	$-0.162$	$3.22\times {10}^{-7}$	$-0.345$	$1.35\times {10}^{-7}$
inf$\_$trim	$0.369$	$1.89\times {10}^{-6}$	$0.696$	$5.09\times {10}^{-7}$
pred$\_$trim	$-0.003$	$5.25\times {10}^{-8}$	$0.017$	$1.49\times {10}^{-8}$
nutri$\_$trim	$0.024$	$5.81\times {10}^{-6}$	$0.115$	$3.71\times {10}^{-5}$
${\sigma }_1$ and ${\sigma }_2$	$1.395$	$4.96\times {10}^{-6}$	1.626	$2.42\times {10}^{-5}$

Equations

$$y_1=X_1{\beta }_1+Z_1{\gamma }_1+{\epsilon }_1,$$

$$y_2=X_2{\beta }_2+Z_2{\gamma }_2+{\epsilon }_2,$$

$$\mathbf{\gamma }=\left(\begin{array}{c}{\gamma }_1\\ {\gamma }_2\end{array}\right)\sim \mathcal{N}\left(0,\mathbf{\Gamma }=\left(\begin{array}{cc}{\Gamma }_1& {\Gamma }_{12}\\ {\Gamma }_{12}^{\top }& {\Gamma }_2\end{array}\right)\right),\mathbf{\epsilon }=\left(\begin{array}{c}{\epsilon }_1\\ {\epsilon }_2\end{array}\right)\sim \left(0,\left(\begin{array}{cc}{\sigma }_1^2{\text{I}}_N& 0\\ 0& {\sigma }_2^2{\text{I}}_N\end{array}\right)\right).$$

$${\Gamma }_1={\sigma }_1^2{\Lambda }_{{\theta }_1}{\Lambda }_{{\theta }_1}^{\top }.$$

$${\Gamma }_2={\sigma }_2^2{\Lambda }_{{\theta }_2}{\Lambda }_{{\theta }_2}^{\top }.$$

$${\gamma }_1={\Lambda }_{{\theta }_1}{u}_1,{\gamma }_2={\Lambda }_{{\theta }_2}{u}_2,$$

$$\mathbf{u}=\left(\begin{array}{c}{u}_1\\ u_2\end{array}\right)\sim \mathcal{N}\left(0,{\Sigma }_{\mathbf{u}}\right),\text{ with }{\Sigma }_{\mathbf{u}}=\left(\begin{array}{cc}{\sigma }_1^2{\text{I}}_{q_1}& {\sigma }_1{\sigma }_2\mathbf{\rho }\\ {\sigma }_1{\sigma }_2{\mathbf{\rho }}^{\top }& {\sigma }_2^2{\text{I}}_{q_2}\end{array}\right).$$

$$\rho =\left(\begin{array}{cc}\text{corr}({\gamma }_1^{\text{I}},{\gamma }_2^{\text{I}})& \text{corr}({\gamma }_1^{\text{I}},{\gamma }_2^{\text{S}})\\ \text{corr}({\gamma }_1^{\text{S}},{\gamma }_2^{\text{I}})& \text{corr}({\gamma }_1^{\text{S}},{\gamma }_2^{\text{S}})\end{array}\right)\text{ and }\mathbf{\rho }=\text{diag}(\rho ,\dots ,\rho ).$$

$$y_1=X_1{\beta }_1+Z_1{\Lambda }_{{\theta }_1}{u}_1+{\epsilon }_1,$$

$$y_2=X_2{\beta }_2+Z_2{\Lambda }_{{\theta }_2}{u}_2+{\epsilon }_2,$$

$$\mathbf{u}=\left(\begin{array}{c}{u}_1\\ u_2\end{array}\right)\sim \mathcal{N}\left(0,{\Sigma }_{\mathbf{u}}\right),\mathbf{\epsilon }=\left(\begin{array}{c}{\epsilon }_1\\ {\epsilon }_2\end{array}\right)\sim \left(0,\left(\begin{array}{cc}{\sigma }_1^2{\text{I}}_N& 0\\ 0& {\sigma }_2^2{\text{I}}_N\end{array}\right)\right).$$

$$\ell (\mathbf{\beta },\mathbf{\theta },\rho ,\mathbf{\sigma }|\mathbf{y})$$

$$\left(\begin{array}{cc}{X}_{\mathbf{\sigma }}^{\top }{X}_{\mathbf{\sigma }}& X_{\mathbf{\sigma }}^{\top }{Z}_{\mathbf{\sigma }\mathbf{\theta }}\\ Z_{\mathbf{\sigma }\mathbf{\theta }}^{\top }{X}_{\sigma }& Z_{\mathbf{\sigma }\mathbf{\theta }}^{\top }{Z}_{\mathbf{\sigma }\mathbf{\theta }}+\sqrt{{\sigma }_1^2{\sigma }_2^2}{\Sigma }_{\mathbf{u}}^{-1}\end{array}\right)\left(\begin{array}{c}{\hat{\mathbf{\beta }}}_{\mathbf{\theta },\rho ,\mathbf{\sigma }}\\ {\mu }_{\mathcal{U}|\mathcal{Y}=\mathbf{y}}\end{array}\right)=\left(\begin{array}{c}{X}_{\mathbf{\sigma }}^{\top }\\ Z_{\mathbf{\sigma }\mathbf{\theta }}^{\top }\end{array}\right)Y_{\mathbf{\sigma }},$$

$$r({\hat{\mathbf{\beta }}}_{\mathbf{\theta },\rho ,\mathbf{\sigma }},{\mu }_{\mathcal{U}|\mathcal{Y}=\mathbf{y}})={\Vert Y_{\mathbf{\sigma }}-X_{\mathbf{\sigma }}{\hat{\mathbf{\beta }}}_{\mathbf{\theta },\rho ,\mathbf{\sigma }}-Z_{\mathbf{\sigma }\mathbf{\theta }}{\mu }_{\mathcal{U}|\mathcal{Y}=\mathbf{y}}\Vert }^2+{\sigma }_1^2{\sigma }_2^2{\mu }_{\mathcal{U}|\mathcal{Y}=\mathbf{y}}^{\top }{\Sigma }_{\mathbf{u}}^{-1}{\mu }_{\mathcal{U}|\mathcal{Y}=\mathbf{y}},$$

$$L_{\mathbf{\theta },\rho ,\mathbf{\sigma }}{L}_{\mathbf{\theta },\rho ,\mathbf{\sigma }}^{\top }=Z_{\mathbf{\sigma }\mathbf{\theta }}^{\top }{Z}_{\mathbf{\sigma }\mathbf{\theta }}+\sqrt{{\sigma }_1^2{\sigma }_2^2}{\Sigma }_{\mathbf{u}}^{-1},$$

$$\left(\begin{array}{cc}{X}_{\mathbf{\sigma }}^{\top }{X}_{\mathbf{\sigma }}& X_{\mathbf{\sigma }}^{\top }{Z}_{\mathbf{\sigma }\mathbf{\theta }}\\ Z_{\mathbf{\sigma }\mathbf{\theta }}^{\top }{X}_{\sigma }& L_{\mathbf{\theta },\rho ,\mathbf{\sigma }}{L}_{\mathbf{\theta },\rho ,\mathbf{\sigma }}^{\top }\end{array}\right)={\left(\begin{array}{cc}{R}_X& 0\\ R_{ZX}& L_{\mathbf{\theta },\rho ,\mathbf{\sigma }}^{\top }\end{array}\right)}^{\top }\left(\begin{array}{cc}{R}_X& 0\\ R_{ZX}& L_{\mathbf{\theta },\rho ,\mathbf{\sigma }}^{\top }\end{array}\right).$$

$$f_{\mathcal{Y}}(\mathbf{y})={\int }_{{\mathbb{R}}^{q_1+q_2}}{f}_{\mathcal{Y},\mathcal{U}}(\mathbf{y},\mathbf{u})d\mathbf{u},$$

$$f_{\mathcal{Y},\mathcal{U}}(\mathbf{y},\mathbf{u})$$

$${\Sigma }_{\mathbf{u}}^{-1}={\tilde{\Sigma }}^{\top }\tilde{\Sigma }.$$

$$\frac{\Vert y_1-X_1{\beta }_1-Z_1{\Lambda }_{{\theta }_1}{u}_1{\Vert }^2}{{\sigma }_1^2}+\frac{\Vert y_2-X_2{\beta }_2-Z_2{\Lambda }_{{\theta }_2}{u}_2{\Vert }^2}{{\sigma }_2^2}+{\mathbf{u}}^{\top }{\Sigma }_{\mathbf{u}}^{-1}\mathbf{u}$$

$$=$$

$$\left(\begin{array}{c}{\hat{\mathbf{\beta }}}_{\mathbf{\theta },\rho ,\mathbf{\sigma }}\\ {\mu }_{\mathcal{U}|\mathcal{Y}=\mathbf{y}}\end{array}\right)=\underset{\mathbf{u},\mathbf{\beta }}{\mathrm{argmin}}g(\mathbf{\beta },\mathbf{u},\mathbf{\theta },\rho ,\mathbf{\sigma })$$

$$Z_{X\Lambda }^{\top }{Z}_{X\Lambda }=\left(\begin{array}{cc}{X}_{\mathbf{\sigma }}^{\top }{X}_{\mathbf{\sigma }}& X_{\mathbf{\sigma }}^{\top }{Z}_{\mathbf{\sigma }\mathbf{\theta }}\\ Z_{\mathbf{\sigma }\mathbf{\theta }}^{\top }{X}_{\sigma }& Z_{\mathbf{\sigma }\mathbf{\theta }}^{\top }{Z}_{\mathbf{\sigma }\mathbf{\theta }}+\sqrt{{\sigma }_1^2{\sigma }_2^2}{\Sigma }_{\mathbf{u}}^{-1}\end{array}\right)\text{ and }{Z}_{X\Lambda }^{\top }{Y}_{\Lambda }=\left(\begin{array}{c}{X}_{\mathbf{\sigma }}^{\top }\\ Z_{\mathbf{\sigma }}^{\top }\end{array}\right)Y_{\mathbf{\sigma }}.$$

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

$$g(\mathbf{\beta },\mathbf{u},\mathbf{\theta },\rho ,\mathbf{\sigma })$$

$$=$$

$$Z_{X\Lambda }^{\top }{Z}_{X\Lambda }={\left(\begin{array}{cc}{R}_X& 0\\ R_{ZX}& L_{\mathbf{\theta },\rho ,\mathbf{\sigma }}^{\top }\end{array}\right)}^{\top }\left(\begin{array}{cc}{R}_X& 0\\ R_{ZX}& L_{\mathbf{\theta },\rho ,\mathbf{\sigma }}^{\top }\end{array}\right),$$

$$L_{\mathbf{\theta },\rho ,\mathbf{\sigma }}{L}_{\mathbf{\theta },\rho ,\mathbf{\sigma }}^{\top }=Z_{\mathbf{\sigma }\mathbf{\theta }}^{\top }{Z}_{\mathbf{\sigma }\mathbf{\theta }}+\sqrt{{\sigma }_1^2{\sigma }_2^2}{\Sigma }_{\mathbf{u}}^{-1}.$$

$$g(\mathbf{\beta },\mathbf{u},\mathbf{\theta },\rho ,\mathbf{\sigma })$$

$$r({\hat{\mathbf{\beta }}}_{\mathbf{\theta },\rho ,\mathbf{\sigma }},{\mu }_{\mathcal{U}|\mathcal{Y}=\mathbf{y}})={\Vert Y_{\mathbf{\sigma }}-X_{\mathbf{\sigma }}{\hat{\mathbf{\beta }}}_{\mathbf{\theta },\rho ,\mathbf{\sigma }}-Z_{\mathbf{\sigma }\mathbf{\theta }}{\mu }_{\mathcal{U}|\mathcal{Y}=\mathbf{y}}\Vert }^2+{\sigma }_1^2{\sigma }_2^2{\mu }_{\mathcal{U}|\mathcal{Y}=\mathbf{y}}^{\top }{\Sigma }_{\mathbf{u}}^{-1}{\mu }_{\mathcal{U}|\mathcal{Y}=\mathbf{y}},$$

$$f_{\mathcal{Y}}(\mathbf{y})$$

Full text

††journal: Journal of LaTeX Templates

1 Introduction

Linear mixed-effects model [Hartley and Rao, 1967, Laird and Ware, 1982, Verbeke, 1997, Hedeker and Gibbons, 2006, Fitzmaurice et al., 2012] has become a popular tool for analyzing univariate multilevel data which arise in many areas (biology, medicine, economy, etc), due to its flexibility to model the correlation contained in these data, and the availability of reliable and efficient software packages for fitting it [Bates et al., 2014, Pinheiro et al., 2007, Littell et al., 1996, Halekoh et al., 2006]. Univariate multilevel data are referred to as observations (or measurements) of a single variable of interest on several levels (school in a village which, in turn, is in a town), while multivariate multilevel data are characterized by multiple variables of interest measured at multiple levels. Examples include exam or test scores recorded for students across time, and multiple items at a single occasion for students in more than one school. Multivariate extension of the (single response variable-based) linear mixed-effects model is, indeed, having increasing popularity as flexible tool for the analysis of multivariate multilevel data [Sammel et al., 1999, Schafer and Yucel, 2002, Wang and Fan, 2010, Jensen et al., 2012].

For the linear mixed-effects model, many methods for obtaining the estimates of the fixed and the random effects have been proposed in the literature. These methods include Henderson’s mixed model equations [Henderson, 1950], approaches proposed by Goldberger [1962] as well as techniques based on two-stage regression, Bayes estimation, etc. For details, see [Searle et al., 1992, Section 7.4c] and Robinson [1991]. Concerning the variance parameters estimation in linear mixed-effects model, the discussed methods in the literature include the ANOVA method for balanced data which uses the expected mean squares approach [Searle, 1995, 1971]. For unbalanced data, Rao [1971] proposed the minimum norm quadratic estimation (MINQUE) method, where the resulting estimates are translation invariant under unbiased quadratic forms of the observations. Lee and Nelder [1998] gave another method of estimating variance parameters using extended quasi-likelihood, i.e. gamma-log generalized linear models. For more details on these parameters’ estimation methods in the linear mixed-effects model, see the paper of Gumedze and Dunne [2011]. Beside all the methods cited earlier, come the Maximum Likelihood (ML) and the Restricted Maximum Likelihood (REML) methods. ML and REML methods are the most popular estimation methods in the linear mixed-effects model [Lindstrom and Bates, 1988]. The main attraction of these methods is that they can handle a much wider class of variance models than simple variance components [Gumedze and Dunne, 2011].

In the multivariate linear mixed-effects model, ML and REML estimates are frequently approached through iterative schemes such as EM algorithm [Meng and Rubin, 1993, Dempster et al., 1977, An et al., 2013, Schafer and Yucel, 2002, Shah et al., 1997]. This avoid the difficulties related to the direct calculating of the parameters’ likelihood, since the random effects are not observed, without ignoring the flexible computationally of these algorithms. Despite the existence of valid theorems which show the asymptotic convergence of the sequences produced by these algorithms toward ML estimates [Dempster et al., 1977], in practice this may not always work exactly as expected.

In this paper, we focus on the multivariate linear mixed-effects model, including all the correlations between the random effects while the marginal residuals are assumed independent homoscedastic with possibly different standard deviation. The class of multivariate mixed-effects models considered here assumes that the random effects and the residuals follow Gaussian distributions, and the dependent variables are continuous. In this model, our approach consists in directly calculating the likelihood of the model’s parameters. This likelihood is used to obtain the ML estimates or the REML estimates through the provided REML criterion. This strategy may explain the high quality of the estimates of both fixed effects parameters and random effects’ variance parameters as well as residual variance parameters. This approach may be viewed as a generalization of the approach proposed by Bates et al. [2014] under the R software [R Core Team, 2015] package named lme4.

2 Multivariate linear mixed-effects model

For the sake of simplicity we focus on the bivariate case ($\displaystyle d=2$) in most of the paper, but the generalization to higher dimensions ($\displaystyle d>2$) is straightforward. Thus, in dimension 2, the model is the following:

[TABLE]

where

[TABLE]

For the sake of simplicity, we write $\displaystyle\boldsymbol{\gamma}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{\Gamma})$ to mean that $\displaystyle\boldsymbol{\gamma}$ is a realization of a random vector which is $\displaystyle\mathcal{N}(\boldsymbol{0},\boldsymbol{\Gamma})$ distributed. For $\displaystyle k\in\{1,2\}$, $\displaystyle\beta_{k}\text{ and }\gamma_{k}$ denote respectively the fixed effects and the random effects vector of covariates, while $\displaystyle\varepsilon_{k}$ is the marginal residual component in the dimension $\displaystyle k$ of the model. $\displaystyle X_{k}$ is a matrix of covariates and $\displaystyle Z_{k}$ a covariates-based matrix of design. dim($\displaystyle X_{k}$)$\displaystyle=N\times p_{k}$ and dim($\displaystyle Z_{k}$)$\displaystyle=N\times q_{k}$, where $\displaystyle N$ is the total number of observations. $\displaystyle p_{k}$ and $\displaystyle q_{k}$ are, respectively, the number of fixed effect related covariates and the number of random effect related covariates in the dimension $\displaystyle k$ of the model. $\displaystyle\boldsymbol{y}=(y_{1}^{\top},y_{2}^{\top})^{\top}$ is the vector of marginal observed response variables of the model. We assume that $\displaystyle\boldsymbol{y}$ is a realization of a random vector $\displaystyle\boldsymbol{\mathcal{Y}}$ and belongs to $\displaystyle\mathbb{R}^{2N}$. The bold symbols represent parameters, or vectors, of multiple dimensions (i.e. $\displaystyle\Gamma_{1}$ concerns dimension $\displaystyle 1$ of the model while $\displaystyle\boldsymbol{\Gamma}$ concerns both dimensions).

$\displaystyle\Gamma_{1}$ and $\displaystyle\Gamma_{2}$ are the variance-covariance matrices of $\displaystyle\gamma_{1}$ and $\displaystyle\gamma_{2}$, respectively. $\displaystyle\Gamma_{1}$ and $\displaystyle\Gamma_{2}$ must be, indeed, positive semidefinite. It is then convenient to express the model in terms of the relative covariance factors, $\displaystyle\Lambda_{\theta_{1}}$ and $\displaystyle\Lambda_{\theta_{2}}$, which are $\displaystyle q_{1}\times q_{1}$ and $\displaystyle q_{2}\times q_{2}$ matrices, respectively. $\displaystyle\Lambda_{\theta_{1}}$ is a block diagonal matrix. Each element in the diagonal of $\displaystyle\Lambda_{\theta_{1}}$ is a lower triangular matrix whose nonzero entries are the components of the vector $\displaystyle\theta_{1}$. That is, $\displaystyle\theta_{1}$ generates the symmetric $\displaystyle q_{1}\times q_{1}$ variance-covariance matrix $\displaystyle\Gamma_{1}$, according to

[TABLE]

Same as $\displaystyle\theta_{2}$ which generates $\displaystyle\Gamma_{2}$ according to

[TABLE]

In Equations 3 and 4, $\displaystyle\sigma_{1}^{2}$ and $\displaystyle\sigma_{2}^{2}$ are the same marginal residual variances used in the model expression (see Equation 2). Using the variance-component parameters, $\displaystyle\theta_{1}$ and $\displaystyle\theta_{2}$, the marginal random effects, $\displaystyle\gamma_{1}$ and $\displaystyle\gamma_{2}$, are expressed as

[TABLE]

such that

[TABLE]

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

[TABLE]

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

[TABLE]

with

[TABLE]

Then the parameters which will be estimated are $\displaystyle\beta_{1},\beta_{2},\sigma_{1}^{2},\sigma_{2}^{2},\theta_{1}$, $\displaystyle\theta_{2}$ and $\displaystyle\rho$.

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.1 ML criterion

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

Theorem 3.1**.**

Suppose that $\displaystyle\boldsymbol{y}=(y_{1}^{\top},y_{2}^{\top})^{\top}$ satisfies the bivariate linear mixed-effects model expressed by Equations (2 and 9), where $\displaystyle\beta_{1},\beta_{2},\sigma_{1}^{2},\sigma_{2}^{2},\theta_{1}$, $\displaystyle\theta_{2}$, $\displaystyle\rho$ are the parameters which need to be estimated, and $\displaystyle\boldsymbol{\beta}=(\beta_{1}^{\top},\beta_{2}^{\top})^{\top}$, $\displaystyle\boldsymbol{\sigma}=(\sigma_{1}^{2},\sigma_{2}^{2})^{\top}$, $\displaystyle\boldsymbol{\theta}=(\theta_{1}^{\top},\theta_{2}^{\top})^{\top}$. Denoting by $\displaystyle Y_{\boldsymbol{\sigma}}=\left(\sqrt{\sigma_{2}^{2}}y_{1}^{\top},\sqrt{\sigma_{1}^{2}}y_{2}^{\top}\right)^{\top}$, $\displaystyle X_{\boldsymbol{\sigma}}=\begin{pmatrix}\sqrt{\sigma_{2}^{2}}X_{1}&\boldsymbol{0}\\ \boldsymbol{0}&\sqrt{\sigma_{1}^{2}}X_{2}\end{pmatrix}$, $\displaystyle Z_{\boldsymbol{\sigma}\boldsymbol{\theta}}=\begin{pmatrix}\sqrt{\sigma_{2}^{2}}Z_{1}\Lambda_{\theta_{1}}&\boldsymbol{0}\\ \boldsymbol{0}&\sqrt{\sigma_{1}^{2}}Z_{2}\Lambda_{\theta_{2}}\end{pmatrix}$, and $\displaystyle\mu_{\boldsymbol{\mathcal{U}}|\boldsymbol{\mathcal{Y}}=\boldsymbol{y}}$ the conditional mean of $\displaystyle\boldsymbol{\mathcal{U}}$ given that $\displaystyle\boldsymbol{\mathcal{Y}}=\boldsymbol{y}$, the log-likelihood of $\displaystyle\boldsymbol{\beta}$, $\displaystyle\boldsymbol{\sigma}$, $\displaystyle\boldsymbol{\theta}$ and $\displaystyle\rho$ given $\displaystyle\boldsymbol{y}$ is expressed as

[TABLE]

where $\displaystyle q=q_{1}+q_{2}$, $\displaystyle\widehat{\boldsymbol{\beta}}_{\boldsymbol{\theta},\rho,\boldsymbol{\sigma}}$ and $\displaystyle\mu_{\boldsymbol{\mathcal{U}}|\boldsymbol{\mathcal{Y}}=\boldsymbol{y}}$ satisfy

[TABLE]

[TABLE]

$\displaystyle L_{\boldsymbol{\theta},\rho,\boldsymbol{\sigma}}$* satisfies*

[TABLE]

and $\displaystyle R_{X}$ satisfies

[TABLE]

Proof.

Denoting by $\displaystyle f_{\mathcal{X}}(.)$ the density function of any random vector $\displaystyle\mathcal{X}$,

[TABLE]

where

[TABLE]

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

[TABLE]

It then comes that $\displaystyle\boldsymbol{u}^{\top}\Sigma_{\boldsymbol{u}}^{-1}\boldsymbol{u}=\|\widetilde{\Sigma}\boldsymbol{u}\|^{2}$ and

[TABLE]

[TABLE]

with

[TABLE]

By setting $\displaystyle p=p_{1}+p_{2}$, $\displaystyle dim(Z_{X\Lambda})=(2N+q)\times(p+q)$ and $\displaystyle S=Im(Z_{X\Lambda})$ is a subspace of $\displaystyle\mathbb{R}^{2N+q}$. $\displaystyle Y_{\Lambda}\in\mathbb{R}^{2N+q}$ and $\displaystyle Z_{X\Lambda}\begin{pmatrix}\widehat{\boldsymbol{\beta}}_{\boldsymbol{\theta},\rho,\boldsymbol{\sigma}}\\ \mu_{\boldsymbol{\mathcal{U}}|\boldsymbol{\mathcal{Y}}=\boldsymbol{y}}\end{pmatrix}$ is the orthogonal projection of $\displaystyle Y_{\Lambda}$ on $\displaystyle S$. Then,

[TABLE]

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

[TABLE]

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

[TABLE]

where

[TABLE]

Thereafter,

[TABLE]

By setting

[TABLE]

and returning to the calculation of $\displaystyle f_{\boldsymbol{\mathcal{Y}}}(\boldsymbol{y})$, it comes

[TABLE]

By setting $\displaystyle\boldsymbol{v}=R_{ZX}(\boldsymbol{\beta}-\widehat{\boldsymbol{\beta}}_{\boldsymbol{\theta},\rho,\boldsymbol{\sigma}})+L_{\boldsymbol{\theta},\rho,\boldsymbol{\sigma}}^{\top}(\boldsymbol{u}-\mu_{\boldsymbol{\mathcal{U}}|\boldsymbol{\mathcal{Y}}=\boldsymbol{y}})$, $\displaystyle d\boldsymbol{u}=\frac{1}{|L_{\boldsymbol{\theta},\rho,\boldsymbol{\sigma}}|}d\boldsymbol{v}$ and

[TABLE]

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

[TABLE]

∎

By profiling out $\displaystyle\boldsymbol{\beta}$, the partially profiled log-likelihood is

[TABLE]

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

[TABLE]

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

Corollary 3.1**.**

Suppose that $\displaystyle\boldsymbol{y}=(y_{1}^{\top},y_{2}^{\top})^{\top}$ satisfies the bivariate linear mixed-effects model expressed by Equations (2 and 9). Taking into account the notations in the Theorem 3.1, the ML estimators $\displaystyle\widehat{\boldsymbol{\beta}}$, $\displaystyle\widehat{\boldsymbol{\sigma}}$, $\displaystyle\widehat{\boldsymbol{\theta}}$, $\displaystyle\widehat{\rho}$ of $\displaystyle\boldsymbol{\beta}$, $\displaystyle\boldsymbol{\sigma}$, $\displaystyle\boldsymbol{\theta}$ and $\displaystyle\rho$ satisfy

[TABLE]

3.2 REML criterion

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

Theorem 3.2**.**

Suppose that $\displaystyle\boldsymbol{y}=(y_{1}^{\top},y_{2}^{\top})^{\top}$ satisfies the bivariate linear mixed-effects model expressed by Equations (2 and 9). Taking into account the notations in the Theorem 3.1, the REML criterion of $\displaystyle\boldsymbol{\sigma}$, $\displaystyle\boldsymbol{\theta}$ and $\displaystyle\rho$ given $\displaystyle\boldsymbol{y}$ is expressed as

[TABLE]

Proof.

[TABLE]

where $\displaystyle\boldsymbol{t}=\boldsymbol{\beta}-\widehat{\boldsymbol{\beta}}_{\boldsymbol{\theta},\rho,\boldsymbol{\sigma}}\implies d\boldsymbol{\beta}=\frac{1}{|R_{X}|}d\boldsymbol{t}$ and $\displaystyle\int_{\mathbb{R}^{p}}(2\pi\sigma_{1}^{2}\sigma_{2}^{2})^{-\frac{p}{2}}\exp\left[-\frac{\left\|\boldsymbol{t}\right\|^{2}}{2\sigma_{1}^{2}\sigma_{2}^{2}}\right]d\boldsymbol{t}=1$. ∎

The REML criterion can also be expressed as

[TABLE]

or as

[TABLE]

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

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:

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

[TABLE]

with

[TABLE]

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 $\displaystyle N$.

We randomly choose $\displaystyle\beta_{1},\beta_{2},\sigma_{1}$ and $\displaystyle\sigma_{2}$ whose values are in the first column of Table 3. $\displaystyle\bar{\boldsymbol{\Gamma}}$ is also randomly chosen such that it is positive definite, with the following form:

[TABLE]

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

[TABLE]

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 , $\displaystyle(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 $\displaystyle 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 1 and Table 2. The Table 1 shows that the asymptotique in the longitudinal data requires not only $\displaystyle n\rightarrow\infty$ and $\displaystyle N\rightarrow\infty$, but also $\displaystyle N/n\rightarrow\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 1, when the total number of observations is $\displaystyle N=1000$, the MSE of $\displaystyle\bar{\boldsymbol{\Gamma}}$ is better for $\displaystyle n=100$, $\displaystyle 0.47$ $\displaystyle(0.03-1.29)$, than for $\displaystyle n=300$, $\displaystyle 0.69$ $\displaystyle(0.06-1.94)$. Observing both Table 1 and Table 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.

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 $\displaystyle(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 $\displaystyle 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 and Table 4. Table 3 contains the empirical means of the estimates with their $\displaystyle 95\%$ CI, and the minimum, the maximum and the average number of iterations. Table 4 contains the empirical relative error of the estimators with their $\displaystyle 95\%$ CI. These results show that in the case of naive initialization, the Cmlme estimators outperform the EM estimators. For example, the component of $\displaystyle\beta_{1}$ which is 14.00 is well estimated by Cmlme, $\displaystyle 14.02$ $\displaystyle(13.27-14.45)$ with an empirical relative error of $\displaystyle 0.02$ $\displaystyle(0.00-0.04)$, but poorly estimated by EM, $\displaystyle-2.05$ $\displaystyle(-4.70--0.40)$ with an empirical relative error of $\displaystyle 1.14$ $\displaystyle(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 whith 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, $\displaystyle\bar{\Gamma}$, is well estimated with advised initializations:

[TABLE]

and

[TABLE]

$\displaystyle\sigma_{\Gamma_{\text{adv}}}^{\text{Cmlme}}$ and $\displaystyle\sigma_{\Gamma_{\text{adv}}}^{\text{EM}}$ contain the standard deviations of the entries of $\displaystyle\bar{\Gamma}_{\text{adv}}^{\text{Cmlme}}$ and $\displaystyle\bar{\Gamma}_{\text{adv}}^{\text{EM}}$, respectively. It seems that the empirical standard deviation of the higher entries of $\displaystyle\bar{\Gamma}$ are bigger with EM than with Cmlme. For example, the standard deviation of $\displaystyle\bar{\Gamma}_{11}=27.77$ is $\displaystyle 7.49$ for EM, but $\displaystyle 5.73$ for Cmlme. Same remark about the standard deviations of $\displaystyle\bar{\Gamma}_{31}$, $\displaystyle\bar{\Gamma}_{32}$ and $\displaystyle\bar{\Gamma}_{33}$, comparing EM and Cmlme. This may be explained by the fact that Cmlme estimators are more consistent than EM estimators. In the case of naive initializations, Cmlme provides a well estimated empirical mean of $\displaystyle\bar{\Gamma}$, when the estimated $\displaystyle\bar{\Gamma}$ provided by EM is very bad (we choose not to show it here).

[TABLE]

$\displaystyle\bar{\Gamma}_{\text{naiv}}^{\text{Cmlme}}$ compared to $\displaystyle\bar{\Gamma}_{\text{adv}}^{\text{Cmlme}}$ and $\displaystyle\sigma_{\Gamma_{\text{naiv}}}^{\text{Cmlme}}$ compared to $\displaystyle\sigma_{\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 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.

5 Application on malaria dataset

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 ($\displaystyle 95\%$) Djènontin et al. [2010] from June $\displaystyle 2007$ to January $\displaystyle 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 ($\displaystyle 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 $\displaystyle 37.5^{\circ}$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 ($\displaystyle>37.5^{\circ}$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 $\displaystyle 3$ months (at $\displaystyle 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 [Anders et al., 2010, McCarthy et al., 2011]. 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 ($\displaystyle\mu$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 5, where Y denotes an antigen which is one of the following:

[TABLE]

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 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 5 are covariates. We then have 8 dependent variables which describe the longitudinal profile (in the child) of the proteins listed in Equation 49.

To illustrate the stability of our approach, we are fitting here a bivariate model to the data, with $\displaystyle\text{IgG1}\_\text{A1}$ and $\displaystyle\text{IgG3}\_\text{A2}$ as dependent variables:

[TABLE]

with

[TABLE]

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 allows to evaluate how much the Cmlme algorithm is sensitive to the starting points. The results are contained in the Table 6.

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

[TABLE]

with an MSE of $\displaystyle 0.0095$.

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 Bates et al. [2014] 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 consitent 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.