EM-REML estimation of covariance parameters in Gaussian mixed models for longitudinal data analysis

This paper presents procedures for implementing the EM algorithm to compute REML estimates of variance covariance components in Gaussian mixed models for longitudinal data analysis. The class of models considered includes random coefficient factors, stationary time processes and measurement errors. The EM algorithm allows separation of the computations pertaining to parameters involved in the random coefficient factors from those pertaining to the time processes and errors. The procedures are illustrated with Pothoff and Roy's data example on growth measurements taken on 11 girls and 16 boys at four ages. Several variants and extensions are discussed.


INTRODUCTION
There has been a great deal of interest in longitudinal data analysis among biometricians over the last decade: see e.g., the comprehensive synthesis of both theoretical and applied aspects given in Diggle et al. [4] textbook. Since the pioneer work of Laird and Ware [13] and of Diggle [3], random effects models [17] have been the cornerstone of statistical analysis used in biometry for this kind of data. In fact, as well illustrated in the quantitative genetics and animal breeding areas, practitioners have for a long time restricted their attention to the most extreme versions of such models viz. to the so called intercept or repeatability model with a constant intra-class correlation, and to the multiple trait approach involving an unspecified variance covariance structure.
Harville [9] first advocated the use of autoregressive random effects to the animal breeding community for analysing lactation records from different parities. These ideas were later used by Wade and Quaas [33] and Wade et al. [34] to estimate correlation among lactation yields produced over different time periods within herds and by Schaeffer and Dekkers [28] to analyse daily milk records.
As well explained in Diggle et al. [3], potentially interesting models must include three sources of variation: (i) between subjects, (ii) between times within a subject and (iii) measurement errors. Covariance parameters of such models are usually estimated by maximum likelihood procedures based on second order algorithms. The objective of this study is to propose EM-REML procedures [1,21] for estimating these parameters especially for those involved in the serial correlation structure (ii).
The paper is organized as follows. Section 2 describes the model structure and Section 3 the EM implementation. A numerical example based on growth measurements will illustrate these procedures in Section 4, and some elements of discussion and conclusion are given in Section 5.

MODEL STRUCTURE
Let y ij be the jth measurement (j = 1, 2, . . . , n i ) recorded on the ith individual i = 1, 2, . . . , I at time t ij . The class of models considered here can be written as follows: where x ij β represents the systematic component expressed as a linear combination of p explanatory variables (row vector x ij ) with unknown linear coefficients (vector β), and ε ij is the random component.
As in [3], ε ij is decomposed as the sum of three elements: The first term represents the additive effect of K random regression factors u ik on covariable information z ijk (usually a (k − 1)th power of time) and which are specific to each ith individual. The second term w i (t ij ) corresponds to the contribution of a stationary Gaussian time process, and the third term e ij is the so-called measurement error. By gathering the n i measurements made on the ith individual such that (1) and (2) can be expressed in matrix notation as and where where G (K×K) is a symmetric positive definite matrix, which may alternatively be represented under its vector form g = vechG. For instance, for a linear regression, g = (g 00 , g 01 , g 11 ) where g 00 refers to the variance of the intercept, g 11 to the variance of the linear regression coefficient and g 01 to their covariance. R i in (5) has the following structure in the general case where σ 2 e I ni = var(e i ), and for stationary Gaussian simple processes, σ 2 is the variance of each w i (t ij ) and H i = {h ij,ij } the (n i × n i ) correlation matrix among them such that h ij,ij = f (ρ, d ij,ij ) can be written as a function f of a real positive number ρ and of the absolute time separation d ij,ij = |t ij − t ij | between measurements j and j made on the individual i.
Classical examples of such functions are the power: f (ρ, d) = ρ d ; the exponential: exp(−d/ρ), and the Gaussian: exp(−d 2 /ρ 2 ), functions. Notice that for equidistant intervals, these functions are equivalent and reduce to a first order autoregressive process (AR1).
R i in (6) can be alternatively expressed in terms of ρ, σ 2 and of the ratio This parameterisation via r = (σ 2 , ρ, λ) allows models to be addressed both with and without measurement error variance (or "nugget" in geostatistics).

EM IMPLEMENTATION
Let γ = (g , r ) be the 3+K(K +1)/2 parameter vector and x = (y , β , u ) be the complete data vector where y = (y 1 , y 2 , . . . , y i , . . . , y I ) and u = (u 1 , u 2 . . . , u i , . . . , u I ) . Following Dempster et al. [1], the EM algorithm proceeds from the log-likelihood L(γ; x) = ln p(x|γ) of x as a function of γ. Here L(γ; x) can be decomposed as the sum of the log-likelihood of u as a function of g and of the log-likelihood of ε * = y − Xβ − Zu as a function of r, L(γ; x) = L(r; ε * ) + L(g; u) + const., Under normality assumptions, the two log-likelihoods in (8) can be expressed as: The E-step consists of evaluating the conditional expectation of the complete data log-likelihood L(γ; x) = ln p(x|γ) given the observed data y with γ set at its current value γ [t] i.e., evaluating the function while the M -step updates γ by maximizing (11) with respect to γ i.e., The formula in (8) allows the separation of Q(γ|γ [t] ) into two components, the first Q u (g|γ [t] ) corresponding to g, and the second Q ε (r|γ [t] ) corresponding to r, i.e., ). (13) We will not consider the maximization of Q u (g|γ [t] ) with respect to g in detail; this is a classical result: see e.g., Henderson [11], Foulley et al. [6] and Quaas [23]. The (k, l) element of G can be expressed as If individuals are not independent (as happens in genetical studies), one has to can be made explicit from (10) as where ) which can be computed from the elements of Henderson's mixed model equations [10,11].
Using the decomposition of R i in (7), this expression reduces to (16) In order to maximize Q ε (r|γ [t] ) in (16) with respect to r, we suggest using the gradient-EM technique [12] i.e., solving the M -step by one iteration of a second order algorithm. Since here E(Ω i ) = σ 2H i , calculations can be made easier using the Fisher information matrix as in [31]. LettingQ = ∂Q/∂r, Q = E(∂ 2 Q/∂r∂r ) the system to solve can be written where ∆r is the increment in r from one iteration to the next.
The expressions forQ andQ are unchanged for models without measurement error; one just has to reduce the dimension by one and use H i in place ofH i .
The minimum of −2L can be easily computed from the general formula given by Meyer [20] and Quaas [23] Here y R #−1 y −θ R #−1 y = [N − r(X)]σ 2 /σ 2 which equals to N − r(X) for σ 2 evaluated at its REML estimate, so that eventually This formula is useful to compute likelihood ratio test statistics for comparing models, as advocated by Foulley and Quaas [5] and Foulley et al. [7,8].

NUMERICAL APPLICATION
The procedures presented here are illustrated with a small data set due to Pothoff and Roy [22]. These data shown in Table I contain facial growth measurements made on 11 girls and 16 boys at four ages (8, 10, 12 and 14 years) with the nine deleted values at age 10 defined in Little and Rubin [14].
The mean structure considered is the one selected by Verbeke and Molenberghs [32] in their detailed analysis of this example and involves an intercept and a linear trend within each sex such that where µ is a general mean, α i is the effect of sex (i = 1, 2 for female and male children respectively), and β i is the slope within sex i of the linear increase with time t measured at age j (t j = 8, 10, 12 and 14 years).
The model was applied using a full rank parameterisation of the fixed effects defined as β = (µ + α 1 , α 2 − α 1 , β 1 , β 2 − β 1 ). Given this mean structure, six models were fitted with different covariance structures. These models are symbolized as follows with their number of parameters indicated within brackets: Distance from the centre of the pituary to the pteryomaxillary fissure (unit 10 −4 m).  Variance covariance structures associated with each of these six models are shown in Table II. Due to the data structure, the power function f (ρ, d) = ρ d (in short POW) reduces here to an autoregressive first order process (AR1) having as correlation parameter ρ 2 .
EM-REML estimates of the parameters of those models were computed via the techniques presented previously. Iterations were stopped when the norm  Tables III and IV. Random coefficient models such as {5} are especially demanding in terms of computing efforts. Models involving time processes and measurement errors require a backtracking procedure [2] at the beginning of the iterative process i.e., one has to compute r [k+1] as the previous value r [k] plus a fraction ω [k+1] of the Fisher scoring increment ∆r [k+1] where r [k] is the parameter vector defined as previously at iteration k. For instance, we used ω = 0.10 up to k = 3 in the case of model 3.
Model comparisons are worthwhile at this stage to discriminate between all the possibilities offered. However, within the likelihood framework, one has to check first whether models conform to nested hypotheses for the likelihood test procedure to be valid.
E.g. model 3 (POW + m-error) can be compared to model 2 (POW), as model 2 is a special case of model 3 for σ 2 e = 0, and also to model 1 (intercept) which corresponds to ρ = 1. The same reasoning applies to the 3-parameter model 4 (intercept + POW) which can be contrasted to model 1 (equivalent to model 4 for ρ = 0) and also to model 2 (equivalent to model 4 for g 00 = 0).
In these two examples, the null hypothesis (H 0 ) can be described as a point hypothesis with parameter values on the boundary of the parameter space which implies some change in the asymptotic distribution of the likelihood ratio statistic under H 0 [29,30]. Actually, in these two cases, the asymptotic null distribution is a mixture 1/2X 2 0 +1/2X 2 1 of the usual chi-square with one degree of freedom X 2 1 and of a Dirac (probability mass of one) at zero (usually noted X 2 0 ) with equal weights. This results in a P-value which is half the standard one i.e., P − value = 1/2Pr[X 2 1 > ∆(−2L) obs ]; see also Robert-Granié et al. [26], page 556, for a similar application.
In all comparisons, model 2 (POW) is rejected while model 1 (intercept) is accepted. This is not surprising as model 2 emphasizes the effect of time separation on the correlation structure too much as compared to the values observed in the unspecified structure (Tab. IV). Although not significantly different from model 1, models 3 (POW + measurement error) and 5 (intercept + linear trend) might also be good choices with a preference to the first one due to the lower number of parameters.
As a matter of fact, as shown in Table III, one can construct several models with the same number of parameters which cannot be compared. There are two models with two parameters (models 1 and 2) and also two with three parameters (models 3 and 4). The same occurs with four parameters although only the random coefficient model was displayed because fitting the alternative model (intercept + POW + measurement error) reduces here to fitting the sub-model 3 (POW + measurement error) due toĝ 00 becoming very small. Incidentally, running SAS Proc MIXED on this alternative model leads tô g 00 = 331.4071,ρ = 0.2395 andσ 2 e = 1.0268 i.e. to fitting model 4 (intercept +   POW). However, since the value of −2L m for model 4 is slightly higher than that for model 3, it is the EM procedure which gives the right answer.

DISCUSSION-CONCLUSION
This study clearly shows that the EM algorithm is a powerful tool for calculating maximum likelihood estimates of dispersion parameters even when the covariance matrix V is not linear in the parameters as postulated in linear mixed models.
The EM algorithm allows separation of the calculations involved in the R matrix parameters (time processes + errors) and those arising in the G matrix parameters (random coefficients), thus making its application to a large class of models very flexible and attractive.
The procedure can also be easily adapted to get ML rather than REML estimates of parameters with very little change in the implementation, involving only an appropriate evaluation of the conditional expectation of u ik u il and of ε * i ε * i along the same lines as given by Foulley et al. [8]. Corresponding results for the numerical example are shown in Tables III and IV suggesting as expected some downward bias for variances of random coefficient models and of time processes.
Several variants of the standard EM procedure are possible such as those based e.g., on conditional maximization [15,18,19] or parameter expansion [16]. In the case of models without "measurement errors", an especially simple ECME procedure consists of calculating ρ [t+1] for σ 2 fixed at σ 2[t] , with σ 2[t+1] being updated by direct maximization of the residual likelihood (without recourse to missing data), i.e., where W i = Z i DZ i + H i with D = G/σ 2 , and which can be evaluated using Henderson's mixed model equations by Finally, random coefficient models can also be accommodated to include heterogeneity of variances both at the temporal and environmental levels [8,24,25,27] which enlarges the range of potentially useful models for longitudinal data analysis.