The heterogeneous linear mixed modelLet Y_i = (Y_i1,…, Y_ini) be the response vector for the n_i measurements of the subject i with i = 1,…, N. The linear mixed model [1] for the response vector Y_i is defined as:

(1)

X_i is a n_i × p design matrix for the p-vector of fixed effects β, and Z_i is a n_i × q design matrix associated to the q-vector of random effects u_i which represents the subject specific regression coefficients. The errors ε_i are assumed to be normally distributed with mean zero and covariance matrix σ²I_ni, and are assumed to be independent from the vector of random effects u_i.

In an homogeneous mixed model [1], u_i is normally distributed with mean μ and covariance matrix D i.e.

(2)

In the heterogeneous mixed model [2–4], u_i is assumed to follow a mixture of G multivariate Gaussians with different means (μ_g)_g=1,G and a common covariance matrix D i.e.

(3)

Each component g of the mixture has a probability π_g and the (π_g)_g=1,G verify the following conditions:

(4)

In this work, we propose a slightly more general formulation of the model described in (1) in which the effect of some covariates may depend on the components of mixture and some of the random effects may have a common mean whatever the component of mixture. Thus, the X_i design matrix is split in X_1i associated with the vector β of fixed effects which are common to all the components and X_2i associated with the vectors γ_g of fixed effects which are specific to the components. The Z_i design matrix is also splitted in Z_1i associated with the vector v_i of random effects following a single Gaussian distribution and Z_2i associated with the vector u_i of random effects following a mixture of Gaussian distributions. The model is then written as:

(5)

where v_i ~ N(0, D_v) and ; given the component g, the conditional distribution of the vector is with .

 

Following the previous works [3,4], we define w_ig the unobserved variable indicating if the subject i belongs to the component g. We have P(w_ig = 1) = π_g. The density for the vector y_i can then be written as:

(6)

Given w_ig, y_i follows a linear mixed model, and the density f(y_i|w_ig = 1) denoted by (_ig is the multivariate Gaussian density with mean E_ig and covariance matrix V_i given by:

(7)

Let now θ be the vector of the m parameters of the model. θ contains ψ with and π the vector of the G &S722; 1 first component probabilities (π_g)_g=1,G&S722;1. Note that π_g is entirely determined by π as . Vec(D) represents the vector of the upper triangular elements of D. The estimates of θ are obtained as the vector that maximizes the observed log-likelihood:

(8)

 

We propose to maximize directly the observed log-likelihood (8) using a modified Marquardt optimization algorithm [9], a Newton-Raphson like algorithm [10]. The diagonal of the Hessian at iteration k, H^(k), is inflated to obtain a positive definite matrix as: with and if i ≠ j. Initial values for λ and η are λ = 0.01 and η = 0.01. They are reduced when H* is positive definite and increased if not. The estimates θ^(k) are then updated to θ(^k+1) using the current modified Hessian H*^(k) and the current gradient of the parameters g(θ^(k)) according to the formula:

(9)

where, if necessary, α is modified to ensure that the log-likelihood is improved at each iteration.

To ensure that the covariance matrix D is positive, we maximize the log-likelihood on the non zero elements of U, the Cholesky factor of D (i.e. U′U = D) [7]. Furthermore, to deal with the constraints on π (4) we use the transformed parameters (γ_g)_g=1,G&S722;1 with:

(10)

Standard errors of the elements of D and (π_g)_g=1,G&S722;1 are computed by the Δ-method [11] while standard errors of the other parameters are directly computed using the inverse of the observed Hessian matrix.

The convergence is reached when the three following convergence criteria are satisfied: , |L^(k)&S722;L^(k&S722;1)| ≤ ε_b and g(θ^(k))′H^(k)&S722;1g(θ^(k)) ≤ ε_d. The default values are ε_a = 10^&S722;5, ε_b = 10^&S722;5 and ε_d = 10^&S722;8.

As the log-likelihood of a mixture model may have several maxima [8], we use a grid of initial values to find the global maximum. The multimodality of the log-likelihood in mixture models has been often discussed and some authors proposed different strategies to choose the set of initial values [12]. However, none of them seems to be optimal in a general way. We have observed, in our experience, that the results were mainly sensitive to initial values of (π_g)_g=1,G&S722;1 and (μ_g)_g=1,G and less sensitive to the other parameters (Vec(U), β and σ) for which estimates of the homogeneous mixed models were good initial values.

A mixture model is estimated with a fixed number of components G, otherwise the number of parameters in the model is unknown. To choose the right number of components, one has to estimate models with different values for G and select the best model according to a test or a criterion. Some works favor a bootstrap approach to approximate the asymptotic distribution of the likelihood ratio test between models with different number of components [13] but this approach is very heavy in particular for mixture models with random effects. Criteria such as Akaike’s Information Criterion (AIC) [14] or Bayesian Information Criterion (BIC) [15] are often preferred. We use these selection criteria to select the optimal number of components.

 

After parameter estimation, mixture models allow to classify subjects according to the G components. The classification is based on the posterior probabilities (π_ig)_g=1,G that the subject i follows each of the G components. Using = (′, ′)′, these probabilities are obtained by the Bayes theorem [2–4] as:

(11)

We then assign to each subject i the component to which he has the highest probability (π_ig)_g=1,G to belong.