Bayesian joint modelling of the mean and covariance structures for normal longitudinal data

Statistics & Operations Research Transactions SORT 31 (2) July-December 2007, 181-200 Statistics & Operations Research Transactions Bayesian joint modelling of the mean and covariance structures for normal longitudinal data c Institut d’Estad´ ıstica de Catalunya [email protected] ISSN: 1696-2281 www.idescat.net/sort Edilberto Cepeda-Cuervo1and Vicente N´ u˜ nez-Ant´ on2 1Universidad Nacional de Colombia, 2Universidad del Pa´ıs Vasco Abstract We consider the joint modelling of the mean and covariance structures for the general antedependence model, estimating their parameters and the innovation variances in a longitudinal data context. We propose a new and computationally efficient classic estimation method based on the Fisher scoring algorithm to obtain the maximum likelihood estimates of the parameters. In addition, we also propose a new and innovative Bayesian methodology based on the Gibbs sampling, properly adapted for longitudinal data analysis, a methodology that considers linear mean structures and unrestricted covariance structures for normal longitudinal data. We illustrate the proposed methodology and study its strengths and weaknesses by analyzing two examples, the race and the cattle data sets. MSC: 62F15; 62J05; 62F10; 62P10 Keywords: Antedependence models; Bayes estimation; Fisher scoring; Gibbs sampling 1 Introduction Continuous longitudinal data consist of repeated measurements on the same subject over time. These measurements are typically correlated and there have been several proposals in the literature to handle stationary or nonstationary correlations and variances, as well as balanced or unbalanced longitudinal data (see, e.g., Laird and Ware, 1982; Diggle et al., 1994 or Zimmerman and N´ u˜ nez-Ant´ on, 2001). Address for correspondence: Vicente N´ u˜ nez-Ant´ on. Departamento de Econometr´ ıa y Estad´ ıstica (E.A.III). Facultad de Ciencias Econ´ omicas y Empresariales. Universidad del Pa´ ıs Vasco-Euskal Herriko Unibertsitatea. Avenida Lehendakari Aguirre 83. E-48015 Bilbao, Spain. Telephone: +34 94 6013749; fax: +34 94 6013754. E-mail: [email protected]. Received: September 2007 Accepted: October 2007 182 Bayesian joint modelling of the mean and covariance structures for normal longitudinal data In the context of the parametric multivariate regression model for longitudinal data and under normality, the response variable for each of the mexperimental units under study, each having nobservations over time, is denoted by Yi=(Yi1,...,Yin), i=1,...,m.Inthisway,thenm ×1 response vector Y=(Y1,...,Ym)contains the responses for all subjects under study, and it is assumed that the Yi’s are independently normally distributed as N(μ μ μ, Σi=σ2In), with μ μ μ=(μ1,...,μ n)=Xβ β βand Inbeing the identity matrix of order n. Here, Xis the n×pdesign matrix containing the set of explanatory variables, and β β βis the set of mean parameters, so that β β β=(β1,...,βp).This model can be formally written as Yi=μ μ μ+  i,with   i∼N(0,Σi=σ2In),(1) As is well known, this model assumes that the errors are independently and normally distributed with mean zero and unknown constant variance σ2. Under the model above, we have that E(Y)=μ μ μ∗=(μ μ μ,...,μ μ μ)and that ΣYis a block-diagonal matrix with diagonal matrix elements given by Σi=σ2In,i=1,...,m. If ij and ik,jk,i=1,...,m, are not independent, then Var(  i)=Σ iis no longer a diagonal matrix and it would be necessary to model and estimate the offdiagonal elements of the covariance matrix. This modelling approach usually requires to impose some constraints on the elements of Σito guarantee its positive definiteness. For example, in stationary Gaussian processes, such as the ones used in Geostatistics, the covariance between two observations is explicitly determined by their correlation function. More specifically, it is modelled as a function of the (Euclidean) distance between these two observations. Moreover, and given that some of the properties of this function are imposed by its spatial structure, only correlation functions belonging to the families where these requirements hold can be considered (see, e.g., Diggle and Verbyla, 1998, or Stein, 1999). Longitudinal data typically consist of several measurements taken over time in each of the experimental units in the sample. It falls into the framework of correlated observations on the same subject and/or experimental unit, and it requires the specification and estimation of both the mean and the covariance structures. Most of the parametric approaches have concentrated on normal linear models (see, e.g., Pourahmadi, 1999 and 2000, or Zimmerman and N´ u˜ nez-Ant´ on, 2001). A central idea to be able to efficiently estimate the covariance matrix was first introduced by Macchiavelli and Arnold (1994) and Macchiavelli and Moser (1997) and it was based on its Cholesky decomposition. This approach has been used for several joint modelling proposals for the mean and covariance structures in the context of longitudinal data (see, e.g., Pourahmadi, 1999 and 2000, Pan and MacKenzie, 2003, 2006 and 2007 or Pan and Ye, 2006). Zimmerman and N´ u˜ nez-Ant´ on (1997) and Zimmerman, N´ u˜ nez-Ant´ on and El Barmi (1998) also proposed a joint modelling of the mean and covariance structures, and N´ u˜ nez-Ant´ on and Zimmerman (2000) addressed the possibility of having random coefficients and other alternative nonstationary models in a joint modelling proposal for Edilberto Cepeda-Cuervo and Vicente N´ u˜ nez-Ant´ on 183 the mean and covariance structures in the context of longitudinal data. In addition, there have been only a few proposals within the Bayesian framework (see, e.g., Cepeda, 2001, Cepeda and Gamerman, 2000 and 2004, Daniels and Pourahmadi, 2002, or Pourahmadi and Daniels, 2002) and all of them proposed specific and restricted parametric structures for the mean, the innovation variances and the autoregressive parameters in the model. Cepeda and Gamerman (2000) proposed a Bayesian methodology for modeling mean and variance heterogeneity, using normal prior distributions for both the mean and variance parameters in the regression model. Cepeda (2001), also using normal prior distributions, extended this methodology to allow for a joint modelling of the mean and covariance structures. These latter results are included in Cepeda and Gamerman (2004). Independently, Daniels and Pourahmadi (2002) and Pourahmadi and Daniels (2002), also proposed the use of normal prior distributions for both the mean and covariance parameters, but they did not include any explicit algorithm to fit the joint mean and covariance model and use a data set (Pourahmadi and Daniels, 2002) and simulations (Daniels and Pourahmadi, 2002) to illustrate their proposals. Moreover, their proposals focused on modelling the covariance structure and did not include simulations or applications where there was a joint modelling approach proposal for the mean and covariance structure. In this paper, we consider the general antedependence model (Gabriel, 1962, Macchiavelli and Arnold, 1994, or Zimmerman and N´ u˜ nez-Ant´ on, 1997), and propose a joint modelling approach for the mean and covariance structures, estimating the mean and autoregressive parameters, and the innovation variances in the longitudinal data context. This general model does not impose any specific or restricted parametric structure on the innovation variances and autoregressive parameters, as was the case in previous proposal. We initially consider a new and computationally efficient classical estimation algorithm based on the Fisher scoring algorithm to obtain the estimators of the parameters. This proposal is very convenient and appealing in many cases, especially in the ones where the number of observational units in the study is large, such as in the examples used here to illustrate our proposals (i.e., the race data and the cattle data). In these cases there is a better agreement between sample regressograms and fitted autoregressive parameters and innovation variances, resulting in a better estimation of the parameters in the mean structure. In addition, we also propose a new and innovative Bayesian methodology based on the Gibbs Sampling (Geman and Geman, 1984), properly adapted for longitudinal data analysis, a methodology that considers linear mean structures and unrestricted covariance structures for normal longitudinal data. This methodology allows the researcher to be able to incorporate relevant prior information in the data analysis, as well as to obtain the parameter estimates when the number of observational units in the study is small. In this specific case, we can also estimate credibility intervals for the parameters of interest in the model. We illustrate the proposed methodology and study its relative strengths and weaknesses when compared to other proposed methods by analyzing two examples, 184 Bayesian joint modelling of the mean and covariance structures for normal longitudinal data the race and the cattle data sets. Moreover and for the cases where no prior information is available to be implemented in our Bayesian methodology proposal, we also include the methodology for the possibility of using noninformative priors. The comparison of the results obtained with the classic methodology proposal and with the noninformative priors Bayesian proposal allows us to be able to evaluate their efficiency. As will be seen in the examples presented here, the estimates obtained under these two alternative proposals are very similar. The paper is organized as follows. Section 2 introduces the general model used in the context of longitudinal data analysis. In Section 3 we introduce the proposed classic methodology for this type of data. Section 4 includes the proposed Bayesian methodology, which is finally applied to the race and cattle data sets in Section 5. Section 6 includes some general conclusions. 2 The General Model As we have already indicated, one of the main issues in the modelling approach we propose requires Var(  i)=Σ i,i=1,...,m, to be nonnegative definite so that its inverse can be efficiently calculated and, in addition, it should also be allowed to have a general form so that it is not too restrictive. Pourahmadi (1999) proposed a general approach where all of these conditions hold. Note that, since observations on different subjects are assumed to be independent and, thus, only within-subject covariance structures need to be considered, we suppress the subscript i(identifying the subjects) when describing these structures. More specifically and following the general model settings presented in Cepeda and Gamerman (2004), let us consider the general antedependence model (Gabriel, 1962 or Zimmerman and N´ u˜ nez-Ant´ on, 1997), where for a given individual having nobservations we have that: Yij −μj= j−1  k=1 φjk(Yik −μk)+νj,ν j∼N(0,σ 2 j),(2) i=1,...,m,j=1,...,n, and that E(Yij)=μj, where μjis assumed to be a linear function of the vector of parameters β β β. In addition, the νj’s are assumed to be mutually independent and, by convention, we set all empty sums to zero, that is 0 k=1zk=0. In this way, (2) can be rewritten in matrix form as ν ν ν=T(Yi−μ μ μ), ν ν ν∼N(0,D),and D=diag(σ2 1,...,σ 2 n),(3) Edilberto Cepeda-Cuervo and Vicente N´ u˜ nez-Ant´ on 185 where ν ν ν=(ν1,...,ν n),μ μ μ=(μ1,...,μ n),T={tij}j=1,...,n i=1,...,n, with tij =⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ 1ifj=i −φij if j<i 0otherwise and Var(ν ν ν)=D=TVar(Yi−μ μ μ)T=TΣiT=TΣT(4) As a direct result of equations (3) and (4), Σcan be indirectly calculated by computing Dand T. In addition, we should point out that the triangular decomposition in equation (4) is unique. Moreover, given that Σis a symmetric matrix if and only if there exists a unique lower triangular matrix T, with ones in the diagonal, and a unique diagonal matrix Dwith positive diagonal entries such that TΣT=D,wealsohavethatΣis positive definite (Pourahmadi, 1999). Therefore, from (3), we have that ˜ Yi=(In−T)˜ Yi+ν ν ν=Φ˜ Yi+ν ν ν, i=1,...,m,(5) where ˜ Yi=(Yi−μ μ μ), and the k-th row of the n×nmatrix Φ=(In−T) contains [n−(k−1)] zeroes and the (k−1) components of the autoregressive parameter vector φ φ φk=(φk1,...,φ k,k−1),k=2,...,n. 3 Classic Methodology Under the assumption that Yi=(Yi1,...,Yin)∼i.i.d.N(μ μ μ, Σ), i=1,...,m, where μ μ μis assumed to depend linearly on β β β(i.e., μ μ μ=Xβ β β,Xbeing the n×pdesign matrix), and Σ−1=TD−1T, the likelihood function is given by L(β β β, Φ,D|Y)∝|D|−m/2exp −1 2(Y−μ μ μ∗)Σ−1 Y(Y−μ μ μ∗), where Y=(Y1,...,Ym)∼N(μ μ μ∗,ΣY)and|Σ|=|T||D||T|=|D|. Note that in the equation above, μ μ μ∗=(E(Y1),...,E(Ym))=(μ μ μ,...,μ μ μ)and ΣYis a block diagonal matrix with diagonal matrix elements given by Σ,sothatΣ−1 Yis a block diagonal matrix with diagonal matrix elements given by Σ−1=TD−1T. Therefore, the log-likelihood function can be written as (β β β, Φ,D|Y)= log L(β β β, Φ,D|Y)∝−mlog |D|−(Y−μ μ μ∗)Σ−1 Y(Y−μ μ μ∗), so that the components of the corresponding score function are given by 186 Bayesian joint modelling of the mean and covariance structures for normal longitudinal data ∂ ∂β β β=−X∗Σ−1 Y(Y−X∗β β β) ∂ ∂φij =−1 2(Y−X∗β β β)∂Σ−1 Y ∂φij (Y−X∗β β β) (i=1,...,n,j=1,...,i−1), where X∗=(X,...,X)is the nm ×pdesign matrix. Thus, we have that Iβ β β,φ =E∂2 ∂β β β∂φ=0 Iβ β β,σ2=E∂2 ∂β β β∂σ2=0 Now, given that the log-likelihood function can be written as =log L∝−mlog |D|− 1 σ2 1 Y∗ 1Y∗ 1−···− 1 σ2 n (Y∗ n−˜ μ μ μn)(Y∗ n−˜ μ μ μn), where Y∗ 1is the m-dimensional vector with i-th component given by (Yi1−μ1), Y∗ k(k=2,...,n)isthem-dimensional vector with i-th component given by (Yik −μk), and ˜ μ μ μk,k=2,...,n,isthem-dimensional vector with i-th component given by φk1(Yi1−μ1)+···+φk,k−1(Yi,k−1−μk−1), we can write ∂ ∂μ1 =1 σ2 1 m  i=1 (Yi1−μ1) ∂ ∂σ2 1 =−m 2σ2 1 +1 2σ4 1 Y∗ 1Y∗ 1 Therefore, the maximum likelihood estimators of μ1and σ2 1are given by ˆμ1=1 mm i=1Yi1 and ˆσ2 1=1 mY∗ 1Y∗ 1. Now, for φ φ φk=(φk1,...,φ k,k−1)andσ2 k,k=2,...,n, and if we let ˜ Xk be an m×(k−1) matrix with columns given by Y∗ 1,...,Y∗ k−1,wehavethat ∂ ∂φ φ φk =1 σ2 1 ˜ X k(Y∗ k−˜ Xkφ φ φk) ∂ ∂σ2 k =−m 2σ2 k +1 2σ4 k (Y∗ k−˜ μ μ μk)(Y∗ k−˜ μ μ μk) Thus, the maximum likelihood estimators of φ φ φkand σ2 k(k=2,...,n) are given by ˆ φ φ φk=(˜ X k˜ Xk)−1(˜ X kY∗ k) ˆσ2 k=1 m(Y∗ k−˜ μ μ μk)(Y∗ k−˜ μ μ μk) (6) The steps of the algorithm used to obtain the maximum likelihood estimators of both the mean and variance parameters follow: Edilberto Cepeda-Cuervo and Vicente N´ u˜ nez-Ant´ on 187 1. Set some arbitrary initial values for φ φ φkand positive initial values for σ2 k,k= 1,...,n. 2. Compute Σand update β β βby using ∂(β β β, Φ,D|Y) ∂β β β=0. 3. Update φ φ φkand σ2 k, by solving the equations ∂(β β β, Φ,D|Y) ∂φ φ φk =0and∂(β β β, Φ,D|Y) ∂σ2 k =0,k=1,...,n. 4. Repeat steps 2 and 3 until convergence. 4 Bayesian Methodology As in the classic approach, we also assume that Yi=(Yi1,...,Yin)i.i.d.∼N(μ μ μ, Σ), i=1,...,m, where μ μ μis assumed to depend linearly on β β β,andΣ−1=TD−1T. Therefore, the likelihood function is given by L(β β β, Φ,D|Y)∝|D|−m/2exp −1 2(Y−μ μ μ∗)Σ−1 Y(Y−μ μ μ∗), where Y=(Y1,...,Ym),|Σ|=|T||D||T|=|D|and Φ=(In−T). If we now let θ θ θ=(β β β, Φ,D), under the Bayesian approach and in order to obtain the posterior distribution for the parameters, we need to assume a prior distribution P(θ θ θ) for θ θ θ. Without loss of generality and for simplicity, we assume independent prior distributions such that β β β∼N(b0,B), φ φ φk∼N(l0k,Σφ φ φk), ψ1=1/σ2 1∼G(α1,λ 1)and ψk=1/σ2 k∼G(αk,λ k)(k=2,...,n), where α α α=(α1,...,α n),λ λ λ=(λ1,...,λ n), Ψ=(ψ1,...,ψ n),andG(r,s), represents the gamma distribution with parameters r>0 and s>0. As usual, another possibility for the prior distribution for θ θ θcould be to assume a noninformative prior distribution such as, for example, assume Jeffreys prior distributions. From Bayes’ theorem, the posterior conditional distribution for β β β,πβ β β=π(β β β|Φ,Ψ,Y), is given by π(β β β|Φ,Ψ,Y)∝exp −1 2(Y−X∗β β β)Σ−1 Y(Y−X∗β β β)−1 2(β β β−b0)B−1(β β β−b0) ∝exp −1 2(β β β−b∗)B∗−1(β β β−b∗),(7) where b∗=B∗(B−1b0+X∗Σ−1 YY)andB∗=(B−1+X∗Σ−1 YX∗)−1. Therefore, we have that πβ β β=π(β β β|Φ,Ψ,Y)∼N(b∗,B∗) and, thus, it would be possible to sample β β βdirectly from πβ β β. That is, values of β β βcan be proposed directly from πβ β βand accepted with probability one. This is the basic description and motivation for the Gibbs sampler (Geman and Geman, 1984). 188 Bayesian joint modelling of the mean and covariance structures for normal longitudinal data Now, given that Σ−1=TD−1T, we have that, for i=1,...,m, ˜ Y iTD−1T˜ Yi=[(In−Φ)˜ Yi]D−1[(In−Φ)˜ Yi] Therefore, by taking into account the independence between individuals and using the equation above, the quadratic form appearing in the log-likelihood function, Q(Y)= (Y−μ μ μ∗)Σ−1 Y(Y−μ μ μ∗), can be rewritten as Q(Y)=1 σ2 1 Y∗ 1Y∗ 1+···+1 σ2 n (Y∗ n−˜ μ μ μn)(Y∗ n−˜ μ μ μn), where ˜ μ μ μk,k=2,...,nis the m-dimensional vector with i-th component given by φk1(Yi1−μ1)+··· +φk,k−1(Yi,k−1−μk−1), Y∗ 1=(Y11 −μ1,Y21 −μ1,...,Ym1−μ1)is the vector of centered observations for all mindividuals at the first time point; Y∗ 2= (Y12−μ2,Y22−μ2,...,Ym2−μ2), is the vector of centered observations for all mindividuals at the second time point; and Y∗ n=(Y1n−μn,Y2n−μn,...,Ymn −μn)is the vector of centered observations for all mindividuals at the n-th time point. In this way, the maximum likelihood function can be written as: L(β β β, Φ,D)∝|D|−m/2exp −1 2σ2 1 Y∗ 1Y∗ 1−···− 1 2σ2 n (Y∗ n−˜ μ μ μn)(Y∗ n−˜ μ μ μn)(8) Thus, if we assume independent normal prior for the φ φ φk’s, we can obtain, from the application of Bayes’ theorem, that the posterior full conditional distribution for φ φ φkis given by π(φ φ φk|β β β, D,Φ−k)∝σ−1 kexp −1 2σ2 k (Y∗ k−˜ μ μ μk)(Y∗ k−˜ μ μ μk)−1 2(φ φ φk−l0k)Σ−1 φ φ φk(φ φ φk−l0k),(9) where Φ−krepresent the parameters in Φexcluding the corresponding ones for φ φ φk. Therefore, the conditional posterior distribution is given by ˜π(φ φ φk|α α α,λ λ λ)∝exp −1 2(φ φ φk−l∗ k)Σ∗−1 k(φ φ φk−l∗ k), where l∗ k=Σ ∗ kΣφ φ φkl0k+1 σ2 k ˜ X k˜ Ykand Σ∗ k=Σ−1 φ φ φk+1 σ2 k ˜ X k˜ Xk−1. From the above, we have that πφ φ φk=π(φ φ φk|β β β, D,Φ−k,Y)∼N(l∗ k,Σ∗ k) (10) So it is possible to sample φ φ φkdirectly from πφ φ φk.Valuesofφ φ φkcan be proposed directly from πφ φ φkand accepted with probability 1. This is the basic description and motivation for the Gibbs sampler (Geman and Geman, 1984). Finally, to be able to sample σ2 k,k=1,2,3,...,n, as we have already seen before, we propose the use of gamma priors for the ψk’s and, thus, from the straight application Edilberto Cepeda-Cuervo and Vicente N´ u˜ nez-Ant´ on 189 of Bayes’ theorem, we obtain gamma posterior distributions for the ψk’s, so that the sampling procedure can be easily handled by using the Gibbs sampler. More concretely, let ˜ Y∗ k=(Y∗ k−˜ μ μ μk) be a random sample of size mfrom a N(0,σ 2 k) distribution (k=1,...,n), with ψk=1/σ2 k.Thatis, ˜ Y∗ krepresents a sample of m individuals at time tk,k=1,...,n. Now, given that the gamma family is closed under sampling, we can assume a gamma prior distribution so that ψk∼G(αk=n0k/2,n0σ2 0k), where n0kis a natural number and σ2 0k>0. Therefore, the posterior distribution of ψk can be directly obtained by using Bayes’ theorem, so that π(ψk|β β β, D,˜ Y∗ k)∝ψ[(n0k+m)/2]−1 kexp{−(n0kσ2 0k+ms2 0k)ψk/2}(11) This expression corresponds to the kernel of the gamma distribution. Thus, we have that π(ψk|β β β(c),φ φ φ(c))=Gn0k+m 2,n0kσ2 0k+ms2 0k 2, where s2 0k=1 mΣm i=1˜ Y∗ k˜ Y∗ k,k=1,...,n. If we decide to assume constant or noninformative priors, the sampling procedure for each of the parameters involved in the estimation process is described below. Given D,Φand a constant prior distribution for β β β, we can sample β β βfrom π(β β β|Φ,Ψ,Y)∝−1 2(β β β−b∗)B∗−1(β β β−b∗),(12) where b∗=B∗(X∗Σ−1 YY)andB∗=(X∗Σ−1 YX∗)−1. Given β β β,Dand a constant prior distribution for the φ φ φk’s, and letting Φ−krepresent the parameters in Φexcluding the corresponding ones for φ φ φk, we can sample the φ φ φk’s (k=1,2,...,n) from (9) by following the procedure below: 1. Sample φ φ φ1from π(φ φ φ1|β β β, D,Φ−1)∝σ−1 1exp −1 2σ2 1 Y∗ 1Y∗ 1). 2. Sample φ φ φ2from π(φ φ φ2|β β β, D,Φ−2)∝σ−1 2exp −1 2σ2 2 (Y∗ 2−˜ μ μ μ2)(Y∗ 2−˜ μ μ μ2). 3. And so on, up to sample φ φ φnfrom π(φ φ φn|β β β, D,Φ−n)∝σ−1 nexp −1 2σ2 n (Y∗ n−˜ μ μ μn)(Y∗ n−˜ μ μ μn). Finally, given D,β β βand a constant prior for the ψk’s, we can sample the ψk’s from π(ψk|β β β(c),φ φ φ(c))=Gm 2,ms2 0k 2, where s2 0kis defined as before. 196 Bayesian joint modelling of the mean and covariance structures for normal longitudinal data covariance structure is more general than theirs. As a way of comparing the estimated cattle weight as a function of time, as in Figure 2, Table 3 shows these estimated values using the estimates obtained from the Bayesian approach. As can be seen, this behaviour is consistent with the behaviour observed in Figure 2. The estimated innovation variances obtained with the proposed Bayesian methodology under noninformative priors with their corresponding estimated standard deviations in parenthesis are ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ˆσ2 1 ˆσ2 2 ˆσ2 3 ˆσ2 4 ˆσ2 5 ˆσ2 6 ˆσ2 7 ˆσ2 8 ˆσ2 9 ˆσ2 10 ˆσ2 11 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠Bayes = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 109.524(30.490) 102.496(42.180) 38.595(15.072) 42.015(13.012) 29.769(8.330) 29.811(8.357) 45.612(12.784) 32.955(9.345) 23.515(6.791) 37.974(11.438) 10.805(3.024) ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ The estimated values of all innovation variances are somehow consistent with the increasing behaviour of variances seen in Table 3 and Figure 2. However, we must be cautious about this issue because these are not the response variances. We will come back to this matter when reporting the estimated variances and correlations for the group A cattle data (see Table 3). If we wish to compare these results with those reported in Pourahmadi (1999) or Cepeda (2001), we could estimate the log-innovation variances. We have done so and estimates are very similar to theirs. Thus, this can be used as a way to indicate that log-innovation variances could be modelled as cubic polynomials (see Pourahmadi, 1999 or Cepeda, 2001). We do not consider it necessary to include these estimated values here. Table 4 shows the estimated autoregressive parameters with the Bayesian approach. The standard deviations obtained with the proposed Bayesian approach show that the autoregressive parameters are different from zero at a 95% credibility level and they are all consistent with the results presented in Pourahmadi (1999). Finally, Table 3 shows the estimated variances and correlations for the group A cattle data obtained with the Bayesian approach. The behaviour of these estimated values is consistent with the one observed in the corresponding sample variances and correlations (not included here for brevity). Edilberto Cepeda-Cuervo and Vicente N´ u˜ nez-Ant´ on 197 Table 4: Estimated autoregressive parameters obtained with the proposed Bayesian methodology for the group A cattle data set. Estimated standard deviations, in parenthesis, are also included. As in the Tmatrix (see Section 2), there are ones on the main diagonal. 1 1.072 1 (0.172) 0.236 0.698 1 (0.110) (0.097) 0.062 −0.416 1.251 1 (0.118) (0.150) (0.181) 0.104 −0.133 −0.015 1.146 1 (0.050) (0.071) (0.059) (0.068) −0.014 −0.200 0.166 0.295 0.789 1 (0.059) (0.087) (0.060) (0.075) (0.039) −0.027 0.017 0.116 −0.278 0.050 1.044 1 (0.085) (0.092) (0.069) (0.099) (0.056) (0.069) −0.182 0.261 0.100 −0.386 0.150 0.538 0.516 1 (0.069) (0.110) (0.059) (0.087) (0.060) (0.053) (0.063) −0.011 0.208 0.166 −0.362 −0.216 −0.204 0.416 1.053 1 (0.077) (0.058) (0.072) (0.062) (0.055) (0.067) (0.037) (0.067) 0.052 −0.306 0.226 −0.147 0.007 0.112 −0.122 0.562 0.617 1 (0.099) (0.103) (0.077) (0.072) (0.058) (0.065) (0.042) (0.083) (0.072) 0.192 −0.330 −0.011 0.269 0.047 −0.292 −0.050 −0.015 0.213 0.907 1 (0.042) (0.051) (0.077) (0.058) (0.048) (0.052) (0.037) (0.057) (0.073) (0.048) 6 Conclusions We have proposed a joint modelling approach for the mean and covariance structures in the context of normal longitudinal data. In the proposals presented here, the mean is modelled in a linear form and the covariance structure is left unrestricted in the sense that no specific parametric model is imposed on it, except for the fact that its modelling makes use of the general antedependence model specification, which has shown to be most useful in practice for nonstationary situations such as the ones present in the data sets analyzed here (see, e.g., Kenward, 1987, Pourahmadi, 1999 and 2000, N´ u˜ nez-Ant´ on and Zimmerman, 2000 or Zimmerman and N´ u˜ nez-Ant´ on, 2001). The proposals include both a classic and a Bayesian approach, allowing for the possibility of having noninformative priors for the latter. The behaviour of the proposed methodology is evaluated by analyzing two data sets and it has proved to be consistent and reasonable when compared to previous and less general proposals. 198 Bayesian joint modelling of the mean and covariance structures for normal longitudinal data Extensions allowing for nonlinear mean structures are being considered at the moment but are beyond the scope of this paper. Acknowledgements Cepeda’s work was supported by a grant from Universidad Nacional de Colombia. N´ u˜ nez-Ant´ on’s work was supported by Ministerio Espa˜ nol de Educaci´ on y Ciencia, FEDER, Universidad del Pa´ ıs Vasco (UPV/EHU) and Departamento of Educaci´ on del Gobierno Vasco (UPV/EHU Econometrics Research Group) under research grants MTM2004-00341, MTM2007-60112 and IT-334-07. The authors thank Dr. Dani Gamerman for helpful comments and suggestions which led to substantial improvement in the presentation of the material in this paper. References Cepeda, E.C. (2001). Variability Modeling in Generalized Linear Models. Unpublished Ph.D. Thesis. Mathematics Institute, Universidade Federal do Rio de Janeiro. Cepeda, E.C. and Gamerman, D. (2000). Bayesian modeling of variance heterogeneity in normal regression models. Brazilian Journal of Probability and Statistics (REBRAPE), 14, 207-221. Cepeda, E.C. and Gamerman, D. (2004). Bayesian modeling of joint regressions for the mean and covariance matrix. Biometrical Journal, 4, 430-440. Daniels, M.J. and Pourahmadi, M. (2002). Bayesian analysis of covariance matrices and dynamic models for longitudinal data. Biometrika, 89, 553-566. Diggle, P.J. and Verbyla, A. (1998). Nonparametric estimation of covariance structure in longitudinal data. Biometrics, 54, 401-415. Diggle, P.J., Liang, K.-Y. and Zeger, S.L. (1994). Analysis of Longitudinal Data. Oxford: Oxford University Press. Gabriel, K.R. (1962). Ante-dependence analysis of an ordered set of variables. Annals of Mathematical Statistics, 33, 201-212. Geman, S. and Geman, D. (1984). Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images. IEEE Transactions on Pattern Analysis and Machine Intelligence, 6, 721-741. Kenward, M.C. (1987). A method for comparing profiles of repeated measurements. Applied Statistics, 36, 296-308. Laird, N.M. and Ware, J.H. (1982). Random effects models for longitudinal data. Biometrics, 38, 963-974. Macchiavelli, R.E. and Arnold, S.F. (1994). Variable order antedependence models. Communications in Statistics. Theory and Methods, 23, 2683-2699. Macchiavelli, R.E. and Moser, E.B. (1997). Analysis of repeated mesurements with ante-dependence covariance models. Biometrical Journal, 39, 339-350. N´ u˜ nez-Ant´ on, V. and Zimmerman, D.L. (2000). Modelling nonstationary longitudinal data. Biometrics, 56, 699-705. Pan, J.X. and MacKenzie, G. (2003). On modelling mean-covariance structures in longitudinal studies. Edilberto Cepeda-Cuervo and Vicente N´ u˜ nez-Ant´ on 199 Biometrika, 90, 239-244. Pan, J.X. and MacKenzie, G. (2006). Regression models for covariance structures in longitudinal studies. Statistical Modelling, 6, 43-57. Pan, J.X. and MacKenzie, G. (2007). Modelling conditional covariance in the linear mixed model. Statistical Modelling, 7, 49-71. Pan, J.X. and Ye, H. (2006). Modelling covariance structures in generalized estimating equations for longitudinal data. Biometrika, 93, 927-941. Pourahmadi, M. (1999). Joint mean-covariance models with applications to longitudinal data: Unconstrained parameterisation. Biometrika, 86, 667-690. Pourahmadi, M. (2000). Maximum likelihood estimation of generalized linear models for multivariate normal covariance matrix. Biometrika, 87, 425-435. Pourahmadi, M. (2002). Graphical diagnostics for modeling unstructured covariance matrices. International Statistical Review, 70, 395-417. Pourahmadi, M. and Daniels, M.J. (2002). Dynamic conditionally linear mixed models for longitudinal data. Biometrics, 58, 225-231. Stein, M.L. (1999). Interpolation of Spatial Data. Some Theory for Kriging. New York: Springer. Wu, W.B. and Pourahmadi, M. (2003). Nonparametric estimation of large covariance matrices of longitudinal data. Biometrika, 90, 831-844. Zimmerman, D.L. (2000). Viewing the correlation structure of longitudinal data through a PRISM. The American Statistician, 54, 310-318. Zimmerman, D.L. and N´ u˜ nez-Ant´ on, V. (1997). Structured antedependence models for longitudinal data. In Modelling Longitudinal and Spatially Correlated Data. Methods, Applications, and Future Directions, T.G. Gregoire, D.R. Brillinger, P.J. Diggle, E. Russek-Cohen, W.G. Warren, and R. Wolfinger (eds.), 63-76. Lecture Notes in Statistics No. 122. New York: Springer-Verlag. Zimmerman, D.L. and N´ u˜ nez-Ant´ on, V. (2001). Parametric modelling of growth curve data: An overview (with comments). Test, 10, 1-73. Zimmerman, D.L., N´ u˜ nez-Ant´ on, V. and El Barmi, H. (1998). Computational aspects of likelihood-based estimation of first-order antedependence models. Journal of Statistical Computation and Simulation, 60, 67-84.