Gaussian estimation of regression and correlation parameters in longitudinal data

Peer review under responsibility of University of Bahrain.

Abstract

The moment estimate of the correlation parameters using generalized estimating equations (GEE) is not guaranteed to exist or to be feasible. We introduce the Gaussian estimation method and show that the estimate of the correlation parameters in longitudinal data setup is asymptotically unbiased and feasible. We derive the large sample properties of the regression and correlation estimates and we illustrate these estimators via real life data example.

Keywords:

• Covariance matrix
• Generalized least squares
• Gaussian estimation
• Generalized estimating equations
• Longitudinal data

1 Introduction

Longitudinal or repeated measures data arise in several areas of research including social, environmental, medical, biological sciences where many experimental units (human, animals, plants, schools, etc.) are measured repeatedly over time. The broad goal of such studies is to assess changes over time within units and relate them to the treatments, interventions and other relevant covariates using the idea of regression. The difficulty in exact specification of the distribution of nonnormal data (binary, count, etc.) precludes the use of maximum likelihood for parameter estimation. However, the idea of quasi-likelihood (QL) for independent observations (Wedderburn, 1974 and Heyde, 1997) requires only a model for the mean and the relationship between the mean and variance. Liang and Zeger (1986) generalized estimating equations (GEE) method extended the framework of QL estimation to dependent longitudinal data by specifying models for the marginal means and variances. They relied on a ‘working’ correlation matrix with few parameters (say, α) to write an analog of the generalized least squares normal equations. Some pitfalls and limitations of the GEE estimator of α were pointed out by Crowder (2001) and other researchers. For example, the moment estimator of the correlation parameter may not exist in some simple cases or it may not be feasible. Despite the list of drawbacks of GEE, researchers insisted on the fact that GEE usually works well in practice even in the misspecified structures (Crowder, 2001). Several articles proposed alternative approaches based on minimizing objective functions instead of solving some ad hoc estimating equations to obtain efficient estimators of β and α (Chaganty, 1997; Chaganty and Shults, 1999; Al-Rawwash, 2001; Crowder, 2001; Al-Rawwash and Pourahmadi, 2006).

Chaganty (1997) relied on partial minimization of the generalized least squares criterion to introduce the quasi-least squares (QLS) estimates where the covariance matrix is a function of both regression and correlation parameters (β, α). The estimating equation for β is precisely the same as GEE, but the estimating equation for α is new and different from all previously developed approaches (see Diggle et al., 1994). Shults and Chaganty (1998) used the QLS method for certain correlation models and showed the importance of the QLS compared to the GEE for unbalanced longitudinal data setup. The fact that the QLS estimator of the correlation parameter is asymptotically biased motivated Chaganty and Shults (1999) to propose a method to eliminate this problem. They proposed using a continuous, one-to-one transformation that depends on the working correlation matrix under the QLS scheme. They presented the method for different correlation structures and obtained a bias corrected estimate of the correlation matrix. Pourahmadi (2000) developed the maximum likelihood estimators for the parameters of a generalized linear model for the covariance matrix and discussed the asymptotic results under normality condition. Al-Rawwash (2005) introduced a nonparametric approach to estimate the parameters of interest as an alternative to the classical parametric techniques. Al-Rawwash and Pourahmadi (2006) introduced the generalized Gaussian likelihood function as their objective function and they used data-driven approach to obtain estimates of the parameters assuming minimal distributional requirements. In this paper, we follow the footsteps of Al-Rawwash and Pourahmadi (2006) and use the generalized Gaussian likelihood function as the objective function to pursue the estimation of the regression and the correlation parameters.

The history of Gaussian estimation (GE) goes back to Whittle (1961) in the time series literature where the parameter estimates are the values that maximize the Gaussian likelihood function for certain correlation structures. Several articles discussed the GE in longitudinal data literature where they focused on obtaining the parameter estimates based on minimizing an objective function even when the observations are not normally distributed (Crowder, 1995; Crowder, 2001; Al-Rawwash, 2001; Al-Rawwash and Pourahmadi, 2006). In this article, we assume only the knowledge of the first two moments and not the exact specification of the distribution of the data. The notion of Gaussian estimation is distinct from GEE in that it minimizes an objective function rather than relying on estimating equations that might produce flawed estimators (Chaganty, 1997; Al-Rawwash, 2001; Crowder, 2001). Several articles discussed the idea of modeling longitudinal data sets including but not limited to Vonesh and Chinchilli (1997), Lipsitz et al. (2000), Wang and Carey (2003), Pourahmadi et al. (2007) and Buzkova and Lamley (2008)".

The organization of the paper is as follows. In Section 2, we introduce the Gaussian estimation and outline a Newton-Raphson iterative method for computing the estimates. Some closed form solutions for the estimates of the correlation parameter α is given in Section 3 when the correlation structure is AR(1). In Section 4, we establish the consistency and asymptotic normality of the regression and the correlation parameter estimates under mild regularity conditions. In Section 5, we provide a real life data analysis of the GE compared to other well known estimation procedures. Finally, we outline the proof of some results in the appendix.

2 The Gaussian estimation

To set the notation, we consider a random vector comprising the measurements taken on a generic subject at times ) and the associated covariates x_j = (x_j1, … ,x_jp)′, j = 1, 2, … ,  n_i, i = 1,2, ⋯ ,m. We plan to make no distributional assumption about Y except for its first two moments. Firstly, the means and variances of the entries of the response are related to the covariates x_j via , var(y_j) = φv(μ_j), where β = (β₁, … , β_p)′ is the regression parameters of primary interest, g is an invertible known link function, v(·) is a known variance function and φ is a dispersion parameter not depending on β. The importance of the parameter φ appears clearly when we handle a nonnormal data set for it expresses the variability beyond the mean. The setup may be carried out for unbalanced longitudinal data, however we focus in this article on the balanced case where n_i = n. Secondly, we use the matrix R = R(α) to model correlations among the measurements on the same subject. The matrix R and its parameters α are usually viewed as nuisance, though in some situations they are of primary interest. Now, if we define μ = μ(β) = (μ₁, … ,μ_n)′ and A(β) = diag(v(μ₁), … ,v(μ_n)), then the covariance matrix of Y is decomposed into

In longitudinal data setup, we assume Y_i, i = 1,2, ⋯ ,m, to be the vector of repeated measures on the ith subject with covariance . The Gaussian estimates (Al-Rawwash and Pourahmadi, 2006) of (β, α) are the minimizers of the objective function(1) where is the vector of partially standardized measurement on the ith subject. Even if G(.,.) is differentiable, its differentiation with respect to β is complicated because of the appearance of β in A_i(β). This problem is avoided by pretending that is a function of β^∗ different from β (sometimes called decoupling). Therefore, if we differentiate (1)(1) with respect to β, α and then substituting β for β^∗ we arrive at:(2) (3) where D(β) = diag(D₁,D₂, ⋯ ,D_m), for 1  ⩽ i ⩽ m, , A(β) = diag(A₁(β),A₂(β), ⋯ , A_m(β)), and ⊗ is the kronecker product. The estimating Eqs. (2)(2) are precisely the GEEs for β, but (3)(3) is different from the estimating equations developed in the literature to estimate the correlation parameter. Solving these two equations for β and α could give the GE of (β,α). Since closed form solution is not available in general, we suggest using the Newton-Raphson iterative algorithm as follows. Firstly, we choose an initial value for β, then we evaluate the following quantities: , , and at . Accordingly, we solve Eq. (3)(3) for , then we compute and for 1 ⩽ i ⩽ m and construct the covariance matrix as well as and . Hence, we update the value of using . Finally, we repeat the process until , then we take as an estimate of β and as an estimate of α. The least squares estimates of β are a good candidate for the initial value of β, while estimates of α must be chosen so that the positive definiteness of the correlation matrix is guaranteed. The estimate of φ depends on the GE results of β and α obtained using Newton-Raphson iterative method which is reduced to The unbiasedness of the GE of (β,α) is not guaranteed when the working correlation structure is misspecified. However, we show in Section 4 that the GE of β and α is consistent when the working correlation matrix is correctly specified.

3 Specific correlation structures

In this section, we present closed form solutions for the correlation parameter of AR(1) correlation structure when the correlation matrix is either misspecified or correctly specified. These can be used in Newton–Raphson algorithm and hence avoid an iterative cycle within this algorithm. The impact of the correlation matrix misspecification on the degrees of estimating equations and the form of their solutions are illuminated. Chaganty (1997) and Chaganty and Shults (1999) considered special correlation structures and they derived closed form solutions for the estimator of the correlation parameter using the the QLS and the C-QLS. The asymptotic bias of the QLS estimator of the correlation parameter is corrected by Chaganty and Shults (1999) using a continuous one-to-one function that transformed the biased correlation matrix to an asymptotically unbiased one. We show that GE of α is asymptotically unbiased and feasible when the number of measurements per subject n is greater than 3 or a slightly modified (conditional) form of the Gaussian estimation is used.

In the following example we assume that the true correlation matrix, denoted by , is unstructured, and we intend to give explicit formulas for the estimate of α by solving (3)(3) for the AR(1) ‘working’ correlation matrices.

Example 1.

Let the working correlation structure of the repeated measurements be the AR(1) matrix R(ρ) = (ρ^∣i−j∣) with the parameter α = ρ∈ (−1, 1). Note thatwhere C₁ is a tridiagonal matrix with 0 on the diagonal and 1 on the first upper and lower diagonals and C₂ = diag(0,1, ⋯ ,1,0). Since the entries of (1 − ρ²)R⁻¹ are quadratic functions of ρ, the estimating Eq. (3)(3) reduces to the quadratic equationwhere , , is the current estimate of β and Z_i is defined in (1)(1) . Its (unique) feasible solution is given byprovided that ∣b_m∣ > ∣a_m∣. This is in sharp contrast to the QLS estimate of the correlation parameter which is known to be feasible all the time (Chaganty, 1997). The differences between the coefficients of the quadratic equations for GE and QLS estimates of α are in the second terms in a_m and b_m, for which with we have . Nevertheless, there are datasets for which ∣b_m∣ < ∣a_m∣ and the quadratic equation have complex roots.

It is important to mention that the working and the true correlation matrices may have the same AR(1) structure which consequently makes corr(Y_i) correctly specified. In this case, the estimating Eq. (3)(3) is cubic in ρ, therefore we shall derive closed formulas for its solutions and study their feasibility.

Example 2

Let both the working and true correlation matrices have AR(1) structures with R(ρ) and R⁻¹(ρ) given in Example 1. Then, (3)(3) reduces to the cubic equation(4) where , and with is defined in (1)(1) . To solve the cubic Eq. (4)(4) , we follow the standard procedure and set whereand for later use define

where the upper sign in γ₁ is used if v is positive and the lower sign if v is negative. The nature and number of feasible roots of (4)(4) depend on the signs of Δ and u as discussed below:

Case 1:	If Δ > 0 and u > 0, then (4)(4) has only one real feasible root, given by
Case 2:	If Δ < 0 and u < 0 , then (4)(4) has three real roots given bywhere , the upper sign in γ is used if v is positive, otherwise, we use the lower sign. Notice that , but the number of feasible roots depends on the data and cannot be decided in advance.
Case 3:	If Δ = 0 and u < 0, then (4)(4) has the unique feasible closed form solution:where the upper sign in is used if v is positive, otherwise we use the lower sign.

3.1 Unique Feasible Solutions

In this section we show that when the number of subjects m is large or when a modified (conditional) Gaussian estimate is used, then a unique feasible correlation parameter estimate satisfying (4)(4) exists.

A number of simulations in Al-Rawwash (2001) corresponding to different values of the triplet (m,n,ρ) revealed that the cubic Eq. (4)(4) has a unique feasible solution for m ⩾ 3. This provided a clue that for m large (4)(4) must have a unique feasible solution in (-1,1), which we show next. From (4)(4) and definitions of a₁ and b₁, it follows that f(−1) < 0 and f(1) > 0. Thus, if we show that f(ρ) is monotone in (-1,1) when m → ∞, then existence of a unique feasible root follows. A sufficient condition for f(ρ) to be monotone is for its derivative f′(ρ) = 3ρ² − 2a₁ρ + b₁ to have at most one real root. But, this happens if and only if its discriminant is nonpositive. Now, applying the law of large numbers (as m → ∞) to a₁ and b₁, it follows thatand the latter is nonpositive if and only if n ⩾ 4.

Another situation that gives rise to a unique feasible correlation parameter estimate for the AR(1) ’working’ correlation structure is obtained by replacing the objective function (1)(1) by that corresponding to the conditional Gaussian likelihood function.

4 Large Sample Properties

This section presents the large sample properties of the Gaussian estimates and . We establish their consistency and asymptotic normality as m → ∞, under certain regularity conditions. Chaganty (1997) showed that the QLS estimate is consistent, but is asymptotically biased and the joint distribution of is asymptotically normal. A corrected version of the QLS to eliminate the bias problem is proposed in Chaganty and Shults (1999) and they considered markov and generalized markov structures to model the unbalanced correlation matrix cases. Thus, compared to the QLS estimate of α, our proposed GE of the correlation parameter is asymptotically unbiased and we intend to compare our results to the QLS as well as the C-QLS estimates. Our proofs are in the spirit of Chaganty (1997); and Sen and Singer (1993, p. 206). Crowder (2001) has presented an alternative proof establishing consistency of the GE, but only asymptotic normality of .

We start by introducing some preliminary facts and notation about the components of (1)(1) . Let Y_i = (y_i1,y_i2, ⋯ , y_in)′,i = 1,2, ⋯ ,m be independent random vectors with E(Y_i) = μ_i(β) and cov(Y_i) = Σ_i(β,α) where the expectation and the covariance are taken under the true probability distribution. Also, let R(α) be the working correlation matrix and let be the true correlation matrix.

Theorem 1

Let θ = (β,α)′, and set(5)

If ∇v_i(θ) and ∇²v_i(θ) exist everywhere and satisfy 0 < E(∇²v_i(θ)) < ∞, then

(a)	E(∇vi(θ)) = 0,
(b)	cov(∇vi(θ)) = ψi(θ),
(c)	E(∇2vi(θ)) = ξi(θ),

wherewith the entries of the above matrices given byAssume that if as m → ∞, we have(6) whereThen,where is a (p ×  1) column vector with one at the j^th row and zero elsewhere, is the true correlation structure and ⊗ is the Kronecker product. It is assumed throughout that components of and are finite.

Next, we present our results on the large sample properties of the GE of (β,α). The regularity conditions (6)(6) needed to establish the consistency and asymptotic normality are similar to those used in the multivariate central limit theorem for independent but not identically distributed random vectors.

Theorem 2

Let θ = (β,α)′ and assume that (6)(6) holds and as δ → 0 we haveThen, the Gaussian estimate is consistent for (β,α) and we conclude thatTheorem 2 establishes the consistency and asymptotic normality of the GE of (β,α) and we outline the proof in the Appendix.

Theorem 3

The marginal distributions of the Gaussian estimates and are asymptotically correlated and

In view of the previous theorems, we notice that if R(α) is correctly specified and is equal to , then we have ψ_i11(θ) = 2ξ_i11(θ), ψ₁₁(θ) = 2ξ₁₁(θ), and the asymptotic covariance of reduces to Also, the asymptotic distribution of depends on the parameters β,α, φ and while the asymptotic distribution of depends on α and but not on β. It also depends on the third and fourth moments of Y_i, which may cause a loss of efficiency if these moments are misspecified. The sandwich-type covariance estimator protects against the impact of misspecification of R(α), however, we may lose some efficiency when using this estimator.

5 Data analysis

Consider the dataset from Potthoff and Roy (1964) that consists of 27 subjects in a dental study (16 boys and 11 girls). The response y_ij is the distance, in mm, from the center of each subjects’ pituitary to pteryomaxillary fissure recorded at the age of 8, 10, 12 and 14. We treat the data set as normally distributed with a mean vector of μ = Xβ and covariance matrix Σ. Also, we consider the models used by Chaganty and Shults (1999) and analyze the data using the GE and compare the results with those in Chaganty and Shults (1999). First, we fit the model:where, x_i1,x_i2 are indicator variables for the two genders: boys and girls, respectively. The covariate x_i3 is the subject’s age at the j^th measurement time. The estimates of (β,α,φ) are computed using the Newton-Raphson iterative method assuming working correlation structures to be:

1.	AR(1) structure,
2.	An unstructured correlation matrix.

Table 1 contains the estimates and the standard errors for the parameters using the GEE, C-QLS and the GE with AR(1) working correlation structure. Similarly, Table 2 contains the estimates when the working correlation matrix is unstructured. The standard error in both Tables were computed using the model-robust, sandwich-type estimator. From Tables 1 and 2, it can be seen that the estimates of the regression parameters are in a reasonable agreement, but the standard errors are slightly different. Overall, the GE and C-QLS have smaller standard errors than the GEE also the GE has a smaller standard errors for all estimates except for the estimate of γ₂ in Table 2. Finally, Tables 3 show the estimates of the working correlation matrices using the GE for the AR(1) and the unstructured correlation matrix.

Table 1 Parameter estimates of the dental data using AR(1) structure.

Parameter	GEE		C-QLS		GE
	Est.	Std.	Est.	Std.	Est.	Std.
β1	16.5946	1.2788	16.5931	1.2781	16.5925	1.2778
β2	17.3213	0.7780	17.3215	0.7776	17.3216	0.7774
γ1	0.7965	0.1050	0.7695	0.1049	0.7696	0.1049
γ2	0.4838	0.0629	0.4837	0.0629	0.4837	0.0629
φ	4.9107		4.9106		4.9105

Table 2 Parameter estimates of the dental data using unstructured correlation.

Parameter	GEE		C-QLS		GE
	Est.	Std.	Est.	Std.	Est.	Std.
β1	16.3236	1.1701	16.0523	1.1288	15.8841	1.1232
β2	17.3973	0.7244	17.4018	0.6972	17.4103	0.6944
γ1	0.7881	0.0983	0.8122	0.0939	0.8276	0.0929
γ2	0.4781	0.0639	0.4770	0.0632	0.4760	0.0642
φ	4.9058		4.9076		4.9101

Table 3 Estimated correlation matrix using GE.

AR(1)				Unstructured
–	0.6082	0.3700	0.2250	–	0.5201	0.6135	0.4873
0.6082	–	0.6082	0.3700	0.5201	–	0.5218	0.7433
0.3700	0.6082	–	0.6082	0.6135	0.5218	–	0.6921
0.2250	0.3700	0.6082	–	0.4873	0.7433	0.6921	–

Notes

Peer review under responsibility of University of Bahrain.

Appendix A

Consider the objective functionwhere, ∣Σ_i∣ = det(Σ_i), θ =  (β,α)′, R(α) is the working correlation matrix, and β is a (p × 1) vector of unknown regression parameters while α is a (q × 1) vector of unknown correlation parameters. In this proof, we use arguments similar to those found in Sen and Singer (1993, p. 206).

Define ε_m(t,θ) by

Under certain regularity conditions and using the Taylor series expansion for around θ and for ∥t∥ ⩽ K, 0 < K < ∞, we can rewrite ε_m(t,θ) as:where θ^∗ is a point on the line joining θ and .

Now, adding and subtracting the term , we get:(A.1) where,Given δ > 0, there exists an integer m₀ = m₀(δ) such thatfor all m ⩾ m₀. Therefore, for sufficiently large m, we have

Hence, using Khintchine strong law of large numbers (Sen and Singer, 1993, p 69) and assuming that E(∇²v_i(θ)) is finite, we havealmost surely as δ→ 0, therefore we have sup {∣η_m(t)∣,∥t∥ ⩽ K} → 0 almost surely as m → ∞.

This will be true since the value of θ^∗ lies on the line joining θ and . NowHence, using the weak law of large numbers, we getTherefore,

In order to find the value of t that minimizes the term ε_m(t, θ), we will find the first derivative of ε_m(t, θ) with respect to t and disregarding the term o_p(1), hence this will give us the following(A.2) Notice that using the definition of ε_m(t,θ), we conclude that will also correspond closely to the Gaussian estimate of θ, which is attained at , therefore we conclude that can be written as follows(A.3) But since and using Theorem 1, (A.2)(A.2) , (A.3)(A.3) and using the weak law of large numbers, we conclude that:(A.4) in probability as m → ∞.

The consistency of Gaussian estimates could be easily concluded from (A.4)(A.4) . Finally, under the assumptions of A.2 ,A.3, A.4 and (6)(6) , the results in Theorem 1, the multivariate central limit theorem and Slutsky’s theorem, we conclude that(A.5) which is the end of the proof.