Research Article Received 22 September 2013,
Accepted 10 April 2014
Published online 8 May 2014 in Wiley Online Library
(wileyonlinelibrary.com) DOI: 10.1002/sim.6194
A note on implementation of decaying product correlation structures for quasi-least squares‡ Justine Shultsa*† and Matthew W. Guerrab This note implements an unstructured decaying product matrix via the quasi-least squares approach for estimation of the correlation parameters in the framework of generalized estimating equations. The structure we consider is fairly general without requiring the large number of parameters that are involved in a fully unstructured matrix. It is straightforward to show that the quasi-least squares estimators of the correlation parameters yield feasible values for the unstructured decaying product structure. Furthermore, subject to conditions that are easily checked, the quasi-least squares estimators are valid for longitudinal Bernoulli data. We demonstrate implementation of the structure in a longitudinal clinical trial with both a continuous and binary outcome variable. Copyright © 2014 John Wiley & Sons, Ltd. Keywords:
decaying product; feasible; generalized estimating equations; longitudinal binary data; quasi-least squares
1. Introduction ( )′ We consider longitudinal studies, in which vectors of outcome measurements Y i = Yi1 , … , Yini and ′ ′ associated covariates xij = (xij1 , … , xijp ) are collected at times Ti = (ti1 , … , tini ) on participant i (i = 1, … , m). The expected value and variance of Yij on subject i can be expressed using a generalized linear model: ( ) ( ) ( ) ( ) (1) E Yij = g−1 xij′ 𝛽 = 𝜇ij and Var Yij = 𝜙h 𝜇ij , respectively. For Bernoulli variables, we assume the logistic model, with g−1 (𝛾) = exp(𝛾)∕(1 + exp(𝛾)) = pij , 𝜙 = 1, and h(𝜇ij ) = pij (1−pij ). For continuous variables, we assume the linear model, with g−1 (𝛾) = 𝛾 and h(𝜇ij ) = 1. We also assume that observations on different subjects are independent, while the repeated measurements are correlated within subjects. This note models the correlation using quasi-least squares (QLS) for unbalanced data [1–3] to implement a decaying product correlation structure. QLS is a two-stage computational approach for estimation of the correlation parameters in the framework of generalized estimating equations (GEE). QLS extends GEE by providing a straightforward approach to implement correlation structures such as the decaying product structure, which previously were not available for GEE. For some structures [1], QLS guarantees an estimated correlation matrix that is feasible (positive-definite), while there is no similar guarantee for the corresponding GEE moment-based estimator. As a result, QLS is sometimes suggested as an alternative approach, should the standard implementation of GEE yield an infeasible estimate of the correlation structure that results in a failure to converge [4]. a Department of Biostatistics and Epidemiology, University of Pennsylvania School of Medicine, Philadelphia, b Division of Biometrics III, Office of Biostatistics, OTS, CDER, FDA, Silver Spring, MD, U.S.A.
PA, U.S.A.
*Correspondence
3398
to: Justine Shults, Department of Biostatistics and Epidemiology, University of Pennsylvania School of Medicine, 423 Guardian Drive, 635 Blockley Hall, Philadelphia, PA 19104, U.S.A. † E-mail: [email protected] ‡ The views expressed in this paper are the authors’ personal views and not necessarily those of the US Food and Drug Administration.
Copyright © 2014 John Wiley & Sons, Ltd.
Statist. Med. 2014, 33 3398–3404
J. SHULTS AND M. W. GUERRA
∏ Decaying product correlations of the form Corr(Yij , Yik ) = k−1 w=j Ciw,w+1 (𝛼) (j < k) are plausible to describe the pattern of association among the repeated measurements on a participant in a longitudinal study. Specifying Cij,j+1 (𝛼) = 𝛼 yields the first order autoregressive or AR(1) structure, for which Corr(Yij , Yik ) = 𝛼 k−j . Setting Cij,j+1 (𝛼) = 𝛼 tij+1 −tij yields the Markov structure, for which Corr(Yij , Yik ) = 𝛼 tik −tij . Setting Cij,j+1 (𝛼) = 𝛼j,j+1 yields the unstructured decaying product, or first-order antedependent ∏k−1 AD(1) structure, for which corr(Yij , Yik ) = w=j 𝛼w,w+1 . For the AD(1) structure 𝛼 = (𝛼1,2 , … , 𝛼N−1,N )′ , where N is the number of pre-planned fixed measurement occasions. [1, 2] implemented the AR(1) and Markov structures for QLS analysis of unbalanced data, while [5] implemented the AR(1), Markov, and AD(1) structures for first-order maximum likelihood (MARK1ML) analysis of longitudinal Bernoulli data. As explained in Section 3.2 of [5], the AR(1) structure is often applied for analysis of longitudinal data that are equally spaced in time, while the more general Markov structure is utilized for unequally spaced measurements. The AD(1) structure can be applied when the number of measurements varies between subjects, but might be most plausible when the subjects have a common set of timings with few intermittent missed visits, that is, most subjects who miss a visit do not return to the study [5]. This note implements the AD(1) structure for QLS. Utilizing the tri-diagonal form of the inverse of the AD(1) structure, it is straightforward to show that the QLS estimators of 𝛼j,j+1 are of the same form as those that were previously derived for the AR(1) structure [1, 2], but restricted to measurement occasions j and j + 1 only. As a result, the estimators are feasible, which means that the QLS estimates always yield a positive definite correlation matrix. Furthermore, subject to conditions that are easily checked [6], the estimators will be valid for longitudinal Bernoulli data, which means that there will be a multivariate distribution with the estimated marginal means and correlation structure. This is important because, while a multivariate (normal) distribution can always be constructed for a given mean and correlation structure for continuous random variables, there is no similar guarantee for Bernoulli variables [7–9]. The AD(1) structure is therefore especially useful for the analysis of longitudinal Bernoulli data because its validity is easily checked and it is a flexible structure with fewer correlation parameters than in a fully unstructured matrix, a feature that is also beneficial for continuous outcomes [10].
2. Implementation of the AD(1) structure for quasi-least squares The following algorithm can be used to implement the decaying product structure for QLS. For simplicity, our presentation assumes a common set of measurement times and that subjects have no intermittent missed visits. However, in the Supporting Materials for the online version of this article (Online Appendix B), we generalize the approach to allow some visits to be missed intermittently. Stage one of QLS: First, select starting values for the stage one estimates 𝛼̃ j−1,j . In our experience, starting values of zero are convenient and perform well, for example, 𝛼̃ j−1,j = 0 ( j = 2, …, N; N = ( )′ ̃ QONE = 𝛼̃ 1,2 , … , 𝛼̃ N−1,N . Next, alternate between Steps A and B: max{ni }). Let 𝜶 Step A: Update the estimate of regression parameter 𝜷 by solving the GEE estimating equation [11] for 𝜷: m ∑
−1∕2 −1 Ri
D′i Ai
) ( ) −1∕2 ( 𝛼̃ QONE Ai Y i − Ui (𝜷) = 0,
(2)
i=1
( )′ ∏ where Ui (𝜷) = 𝜇i1 , 𝜇i2 , · · · , 𝜇ini , Di = 𝜕Ui (𝜷)i ∕𝜕𝜷, Ri [j, k] = k−1 w=j 𝛼w,w+1 (j < k), Ri [k, j] = Ri [j, k], )) ( ( and Ai = diag h(𝜇i1 ), · · · , h 𝜇inii . Step B: Update the estimate of the correlation by obtaining a solution (for 𝜶) to the stage one estimating ̂ for 𝜷: equation evaluated at the current estimate 𝜷 ( ) 𝜕R−1 (𝜶) ( ) i ̂ ̂ = 0 (j = 2, … , N), Zi′ 𝜷 Zi 𝜷 𝜕𝛼 j−1,j i=1
m ∑
where
(
) ( )′ Y i − Ui (𝜷) = zi1 , … , zini
3399
−1∕2
Zi (𝜷) = Ai
(3)
is the vector of Pearson residuals on subject i. Copyright © 2014 John Wiley & Sons, Ltd.
Statist. Med. 2014, 33 3398–3404
J. SHULTS AND M. W. GUERRA
Stage two of QLS: After convergence in the parameters in stage one, we next obtain a consistent estimate of 𝜶 by solving the stage two estimating equation for 𝜶: m ∑
( trace
i=1
𝜕R−1 (𝜶) || i Ri (𝜶) | 𝜕𝛼j−1,j ||𝜶=𝛼̃ QONE
) = 0 (j = 2, … , N).
(4)
( )′ ̂1,2 , … , 𝛼 ̂N−1,N . We then solve the GEE estimating We refer to the solution of (4) as 𝛼 ̂QLS = 𝛼 ( ) equation (2) at R−1 𝛼 ̂QLS , in order to obtain the final estimate of 𝜷. i To simplify (3) and (4), we utilize the simple tri-diagonal form for the inverse of Ri (𝜶) of ) the inverse of for an ni × ni (ni ∈ {2, … , N}) decaying product structure. The elements ( −1 2 [j, j] = ni × ni decaying product AD(1) matrix Ri (𝜶) are [12]: Ri [1, 1] = 1∕ 1 − 𝛼1,2 ; R−1 i ( ) (( )( )) 2 2 2 2 −1 1 − 𝛼j−1,j 𝛼j,j+1 ∕ 1 − 𝛼j−1,j 1 − 𝛼j,j+1 (j = 2, … , ni − 1 for ni > 2); Ri [j, j + 1] = ) ) ( ( 2 2 −𝛼j,j+1 ∕ 1 − 𝛼j,j+1 (j = 1, … , ni − 1); R−1 ; and R−1 [n , n ] = 1∕ 1 − 𝛼 [j, j′ ] = 0 for |j − j′ | > 1. i i i i ni −1,ni The derivative of the inverse of ni × ni (ni ⩾ 2) decaying product matrix Ri (𝜶) with respect to 𝛼j−1,j is then easily shown to equal an ni × ni matrix with the following 2 × 2 sub-matrix in rows j − 1 to j and columns j − 1 to j: 𝜕R−1 (𝜶) [ i 𝜕𝛼j−1,j
) ( 2 ⎛ ⎞ − 1 + 𝛼 2𝛼 j−1,j 1 j−1,j ⎟ ⎜ ( ) j − 1 ∶ j, j − 1 ∶ j = ( )2 ⎜ 2 ⎟ 2𝛼j−1,j 2 ⎝ − 1 + 𝛼j−1,j ⎠ 1 − 𝛼j−1,j ]
(j = 2, … , ni )
(5)
and, for ni > 2, all other elements equal to 0. The solutions to the stage one and stage two estimating equations (3) and (4) can then be directly obtained as follows: √ ) m ( m ( m ( )2 ∑ )2 ∑ ∑ ̂z2ij + ̂z2ij−1 − ̂zij + ̂zij−1 ̂zij − ̂zij−1 𝛼 ̃j−1,j =
i=1
i=1
∑ m
2
i=1
( j = 2, … , N)
(6)
̂ziĵzij−1
i=1
and 𝛼 ̂j−1,j =
2𝛼̃ j−1,j 2 1 + 𝛼̃ j−1,j
( j = 2, … , N).
(7)
The summations in (6) are taken over subjects who have measurements at both visits j and j - 1, so that 𝛼 ̂1,2 (for example) is based on ̂zi1 and ̂zi2 from subjects who were present at visits 1 and 2. The estimators given in (6) and (7) are of the same form as the stage one and two estimators for the AR(1) structure [1, 2], respectively, for a subject with measurements at visits j − 1 and j only. Therefore, the proofs that were provided earlier [1, 2] are applicable to show that both −1 < 𝛼 ̃j−1,j < 1 and −1 < 𝛼 ̂j−1,j < 1 ∀j = 2, … , N; as a result, the estimated AD(1) matrix is positive definite, so that the QLS estimators are feasible for this structure. In addition, a valid multivariate distribution can be constructed ̂ and 𝜶 ̂ QLS as long as the estimates satisfy the Prentice using model (6) in [5] for the given QLS estimates 𝜷 constraints [6] for the consecutive correlations: ̂j−1,j ⩽ Upperi ( j − 1, j), Loweri ( j − 1, j) ⩽ 𝛼
(i = 1, … , m ; j = 2, … , ni )
(8)
3400
{ ( )1∕2 ( )−1∕2 } where Loweri (j − 1, j) = max − w , Upperi (j − 1, j) = min ̂ ij−1 w ̂ ij ,− w ̂ ij−1 w ̂ ij {( )1∕2 ( )−1∕2 } ( )−1 wij , w ̂ ij−1 ∕̂ wij pij 1 − ̂ w ̂ ij−1 ∕̂ pij . In general, the QLS estimators for all , and w ̂ ij = ̂ decaying product structures (Markov, AR(1), AD(1)) will be valid for Bernoulli data as long as the estimated consecutive correlations satisfy the standard Prentice constraints (8) for a bivariate Bernoulli distribution [5]. For other structures, for example, the exchangeable structure with Corr(Yij , Yik ) = 𝛼 for all i, j, k, satisfaction of the constraints (8) must be checked for all pairs. Copyright © 2014 John Wiley & Sons, Ltd.
Statist. Med. 2014, 33 3398–3404
J. SHULTS AND M. W. GUERRA
2.1. Demonstration We next implement the AD(1) structure for a QLS analysis of a longitudinal trial to compare Effexor with lithium in the treatment of bipolar type II depression. Approximately 30% of measurements were randomly dropped in an earlier demonstration, but here we consider the full data set from the Effexor study that was described previously [3, Section 9.3]. We implement model (1) with logit and identity links for a binary and continuous outcome, respectively. In addition, xij′ 𝛽 = 𝛽0 + 𝛽1 I(Effexor) + 𝛽2 visit + 𝛽3 I(Effexor) × visit
(9)
where I(Effexor) takes value 1 or 0 for the 26 patients treated with Effexor or lithium, respectively, and visit = 1, … , ni (ni ∈ {2, … , 8}). There were no subjects with intermittent missing measurements and the mean (median) number of measurements was 6.4 (8.0). If the regression parameter 𝛽3 that is associated with the Effexor by visit interaction differs significantly from zero, this will indicate that the mean change over time in HAMD-17 depression scores (continuous outcome) or log odds of depression (Bernoulli outcome) differs significantly between the two treatment groups. Depression was defined as having a HAMD-17 score > 7. Table I displays the results of the analysis for the HAMD-17 scores, while Table II shows the results for the binary outcome. Several working correlation structures were implemented in the analysis. Results for the fully unstructured matrix are not shown because convergence was not achieved for this structure, for either the binary or the continuous outcome. The DBAR, CONE, and CTWO Rotnizky-Jewell criteria [13] that were assessed in [14] were used to aid in the selection of a working correlation structure. Values of (DBAR, CONE, CTWO) that are closer to (0, 1, 1) indicate better fit. For the continuous outcome, the DBAR criterion selected the AD(1) structure, while the CONE and CTWO criteria selected the AR(1) structure. However, the fit was similar and good for both the AD(1) and AR(1) structures, worse for the exchangeable structure, and very poor for the identity structure. For
Table I. Estimated values of the regression and correlation parameters in the analysis of data from a clinical trial to compare Effexor versus lithium with respect to mean change in HAMD-17 depression scores. Structure AR(1) EXC AD(1) IND
𝛽̂0
𝛽̂1
𝛽̂2
𝛽̂3
𝛼 ̂QLS
pval
DBAR
CONE
CTWO
22.57 21.99 21.93 21.71
−0.33 −1.04 −0.53 −0.65
−1.21 −1.30 −1.06 −1.17
−0.69 −0.54 −0.67 −0.67
0.76 0.62 ** 0
0.05 0.11 0.07 0.11
0.31 0.95 0.29 10.19
0.95 1.50 0.99 3.14
1.21 2.95 1.28 15.46
Results are shown for the AR(1), exchangeable (EXC), AD(1), and identity structures. The pvalue (pval) for the two-sided test of Ho ∶ 𝛽3 = 0 and the CONE, CTWO, and DBAR criteria for selection of a working correlation are also provided. ̂1,2 = 0.68, 𝛼 ̂2,3 = 0.81, 𝛼 ̂3,4 = 0.62, 𝛼 ̂4,5 = 0.69, 𝛼 ̂5,6 = 0.82, 𝛼 ̂6,7 = 0.80, 𝛼 ̂7,8 = 0.90 ** 𝛼
Table II. Estimated values of the regression and correlation parameters in the analysis of data from a clinical trial to compare Effexor versus lithium with respect to change in the log odds of depression. Structure AR(1) EXC AD(1) IND
𝛽̂0
𝛽̂1
𝛽̂2
𝛽̂3
𝛼 ̂QLS
Lowers
Upper
pval
DBAR
CONE
CTWO
0.47 0.50 0.48 0.43
0.42 −0.28 0.55 0.14
−0.22 −0.28 −0.23 −0.22
−0.33 −0.31 −0.36 −0.27
0.55 0.41 * 0
−0.04 −0.02 ** ***
0.76 0.13 ** ***
0.10 0.30 0.07 0.18
0.42 5.69 0.53 5.02
1.30 2.42 1.37 2.74
2.02 9.52 2.27 9.51
Copyright © 2014 John Wiley & Sons, Ltd.
Statist. Med. 2014, 33 3398–3404
3401
Results are shown for the AR(1), exchangeable (EXC), AD(1), and identity structures. The p-value (pval) for the two-sided test of Ho ∶ 𝛽3 = 0 and the CONE, CTWO, and DBAR criteria for selection of a working correlation are ̂QLS are also provided. also provided. The lower (Lower) and upper (Upper) Prentice constraints (8) for 𝛼 ̂1,2 = 0.60(−0.64, 0.74); 𝛼 ̂2,3 = 0.68(−0.36, 0.74); 𝛼 ̂3,4 = 0.54 (-0.2,0.74); 𝛼 ̂4,5 = 0.54(−0.20, 0.74); 𝛼 ̂5,6 = *𝛼 0.29(−0.11, 0.74); 𝛼 ̂6,7 = 0.51(−0.06, 0.74); 𝛼 ̂7,8 = 0.59(−0.03, 0.74) ** The lower and upper Prentice constraints (8) are provided in parentheses following each aforementioned estimate. *** The Prentice constraints (8) are always satisfied for the identity structure.
J. SHULTS AND M. W. GUERRA
example, DBAR was 0.29 and 0.31 for the AD(1) and AR(1) structures, versus 0.95 and 10.19 for the exchangeable and identity structures, respectively. Assessment of the sample variogram [15, p. 49] also suggested that there is a decay in the intra-subject correlations that is not compatible with an identity or an exchangeable working structure. The sample variogram is described in the Supporting Materials for the online version of this article (Online Appendix A). For the binary outcome, all criteria selected the AR(1) structure. However, the fit was similar for the AR(1) and AD(1) structures and was much worse for both the exchangeable and identity structures. For example, DBAR was 0.42 and 0.53 for the AR(1) and AD(1) structures versus 5.69 and 5.02 for the exchangeable and identity structures, respectively. It is straightforward to show [9, Table 1] that the upper Prentice constraint for model (9) for both the AR(1) and AD(1) structures is { ( ) ( )} min exp −|𝛽̂2 |∕2 , exp −|𝛽̂2 + 𝛽̂3 |∕2 , while the lower Prentice constraint is always negative. The Prentice constraints were not violated for the AR(1) and AD(1) structures, but the estimate 𝛼 ̂QLS = 0.41 clearly violated the estimated Prentice constraints (−0.02, 0.13) for the exchangeable structure. The exchangeable structure should therefore be ruled out according to the rule out criterion [14]. The Prentice constraints will never be violated for the identity structure because the lower bound for the constraint is always negative, while the upper bound is positive. However, the poor fit of the identity structure suggests that it should also be ruled out in favor of the AR(1) or AD(1) structures, which had better (and similar) fit. The results of the analysis were similar for all structures, although the results for the AR(1) and AD(1) structures were most suggestive of a beneficial impact of the intervention. For example, if the criterion for statistical significance was a p-value ⩽ 0.1, then the results for both the continuous and binary outcomes would be significant only for the AD(1) and AR(1) structures. However, it is important to acknowledge that we implemented sandwich-type standard errors that can be underestimated in small samples [16]. Work is needed to assess small sample corrections to sandwich-type estimators of covariance for QLS, similar to assessments made for GEE [17]. The results in Tables 1 and 2 were based on a full data set with no intermittent missed measurements. However, in the Supporting Materials for the online version of this article (Online Appendix B), we present the results for the same data set, but with some randomly dropped measurements and intermittently missed visits. The Supporting Materials for the online version of this article (Online Appendix C) present the results of a small simulation study to assess the implementation of the AD(1) structure for QLS when model (1) is correctly specified and several working structures (Markov, exchangeable, unstructured, and AD(1)) are implemented in the analysis. We considered settings for which the Markov structure was the only true structure, and for which the AD(1) and unstructured matrices were the only true structures. We also assessed the conditions under which violation of the Prentice constraints (8) is more likely to occur. The simulations demonstrated that correctly specifying the working correlation structure, whether Markov or AD(1), resulted in greater small sample efficiency in estimation of the regression parameters. In addition, failure to converge was not an issue when fitting models to the longitudinal Bernoulli data that we simulated using the algorithms in [18, 19]. The simulations also suggested that violation of the Prentice constraints is not likely for a correctly specified working structure, unless the true correlation is very close to a boundary value for the constraints; even then, the violation is not likely to be severe (as assessed by the distance between the estimated correlation parameter and the estimated closest boundary value), especially for increasingly larger sample sizes. However, if the correlation structure is incorrectly specified, then there can be a severe violation of constraints that is not alleviated by increasing sample size. For example, we considered settings for which 𝛼̂ for the incorrectly specified exchangeable structure converges in probability to a value that exceeds the upper boundary value of the Prentice constraints; for these scenarios, the percentage of simulation runs that resulted in a violation of constraints for the incorrectly specified exchangeable structure increased to 100% when the group size increased to m = 60.
3. Conclusion
3402
This note implemented the AD(1) structure for QLS and demonstrated the results in analysis of continuous depression scores and Bernoulli depression status in a randomized clinical trial. Expanding the list of correlation structures available for GEE and QLS allows for more careful modeling of the correlation structure for GEE, which is helpful in avoiding the serious loss in efficiency that can result from Copyright © 2014 John Wiley & Sons, Ltd.
Statist. Med. 2014, 33 3398–3404
J. SHULTS AND M. W. GUERRA
incorrect application of an identity structure [20–23]. In addition, unlike the estimators of the regression parameters, the GEE and QLS estimators of the association parameters may fail to be consistent when the true correlation structure is misspecified; in this situation, Sutradhar and Das [21] suggest that the sandwich-based covariance matrix should be evaluated at the limiting values of the association parameters. As a result, application of the sandwich covariance matrix to correct the standard errors for potential misspecification of the true structure can also require careful assessment of the correlation structure. The AD(1) structure is more flexible than the AR(1) or exchangeable structures that assume a constant value for the consecutive and all pairwise correlations, respectively. The Markov structure that depends on the spacing of measurements might be more appropriate when the temporal spacing of measurements is highly irregular, either due to study design (as in [24]) or to many intermittently missing visits. The assumption of a decaying product pattern in the correlations is not as flexible as the fully unstructured matrix, which can be implemented for balanced data (all ni = n) for QLS [3, Section 5.6]. However, the increased flexibility of the fully unstructured matrix can come at a cost of problems with respect to convergence, which was the case for both outcomes of this study, when implementation of the unstructured matrix was attempted with GEE. The analysis was conducted using the xtqls command [4] in Stata 13.0 that we updated to fit the AD(1) structure. We also updated the qls function utilized in the qlspack package [25] to allow for implementation of the AD(1) structure in R. Stata code to recreate the analysis results in Section 2.1 and Online Appendices A and B is provided in the Supporting Materials for the online version of the article (Online Appendix D). R code to recreate the analysis results in Section 2.1 is provided in Online Appendix E.
Acknowledgements The authors are very grateful to the associate editor and two anonymous reviewers for their careful review and thoughtful suggestions, which allowed us to make considerable improvements to the manuscript.
References
Copyright © 2014 John Wiley & Sons, Ltd.
Statist. Med. 2014, 33 3398–3404
3403
1. Shults J, Chaganty NR. Analysis of serially correlated data using quasi-least squares. Biometrics 1998; 54(4): 1622–1630. 2. Chaganty NR, Shults J. On eliminating the asymptotic bias in the quasi-least squares estimate of the correlation parameter. Journal of Statistical Planning and Inference 1999; 76:145–161. 3. Shults J, Hilbe J. Quasi-Least Squares Regression. Chapman & Hall/CRC Press Monographs on Statistics & Applied Probability: Boca Raton, 2014. 4. Shults J, Ratcliffe S, Leonard M. Improved generalized estimating equation analysis via xtqls for implementation of quasileast squares in Stata. Stata Journal 2007; 7:147–166. 5. Guerra MW, Shults J, Amsterdam J, Ten-Have T. The analysis of binary longitudinal data with time-dependent covariates. Statistics in Medicine 2012; 10:931–948. 6. Prentice RL. Correlated binary regression with covariates specific to each binary observation. Biometrics 1988; 44: 1033–1048. 7. Ziegler A, Vens M. Generalized estimating equations: notes on the choice of the working correlation matrix. Methods of Information in Medicine 2010; 5:421–425. 8. Molenberghs G. Editorial. “Generalized estimating equations: notes on the choice of the working correlation matrix”. Methods of Information in Medicine 2010; 49:419–420. 9. Shults J. Discussion of “Generalized estimating equations: notes on the choice of the working correlation matrix” continued. Letter to the editor. Methods of Information in Medicine 2011; 50:96–99. 10. Westgate PM. A bias correction for covariance estimators to improve inference with generalized estimating equations that use an unstructured correlation matrix. Statistics in Medicine 2012; 35(13):2850–2858. 11. Liang KY, Zeger SL. Longitudinal data analysis using generalized linear models. Biometrika 1986; 73:13–22. 12. Zimmerman DL, Nu˜nez-Antón VA. Antedependence Models for Longitudinal Data. Chapman and Hall/CRC Press: Boca Raton, 2012. 13. Rotnitzky A, Jewell NP. Hypothesis testing of regression parameters in semiparametric generalized linear models for cluster correlated data. Biometrika 1990; 77:485–497. 14. Shults J, Sun W, Tu X, Kim H, Amsterdam J, Hilbe JM, Ten-Have T. A comparison of several approaches for choosing between working correlation structures in generalized estimating equation analysis of longitudinal binary data. Statistics in Medicine 2009; 28:2338–2355. 15. Diggle P, Heagerty P, Liang KY, Zeger S. Analysis of Longitudinal Data. Oxford University Press: New York, 2002. 16. Mancl LA, DeRouen TA. A covariance estimator for GEE with improved small-sample properties. Biometrics 2001; 57(1):126–134. 17. Teerenstra S, Bing L, Preisser JS, Achterberg T, Borm G. Sample size considerations for GEE analyses of three-level cluster randomized trials. Biometrics 2010; 66(4):1230–1237.
J. SHULTS AND M. W. GUERRA 18. Guerra MW, Shults J. A note on the simulation of overdispersed random variables with specified marginal means and product correlations. The American Statistician. In Press; Available from: http://www.tandfonline.com/doi/ [Accessed on 18 Febuary 2014]. 19. Preisser JS, Qaqish BF. A comparison of methods for simulating correlated binary variables with specified marginal means and correlations. Journal of Statistical Computation and Simulation. In Press; Available from: http://www.tandfonline. com/ [Accessed on 18 February 2014]. DOI: abs/10.1080/00949655.2013.818148. 20. Fitzmaurice G. A caveat concerning independence estimating equations with multivariate binary data. Biometrics 1995; 51(1):309–317. 21. Sutradhar B, Das K. On the efficiency of regression estimators in generalised linear models for longitudinal data. Biometrika 1999; 86(2):459–465. 22. Sutradhar B, Das K. On the accuracy of efficiency of estimating equation approach. Biometrics 2000; 56(2):622–625. 23. Wang Y, Carey V. Working correlation structure misspecification, estimation and covariate design: implications for generalised estimating equations performance. Biometrika 2003; 90(1):29–41. 24. Preisser J, Arcury T, Quandt S. Detecting patterns of occupational illness clustering with alternating logistic regressions applied to longitudinal data. American Journal of Epidemiology 2003; 158(5):495–501. 25. Xie J, Shults J. Implementation of quasi-least squares With the R package qlspack. UPenn Biostatistics Working Papers. Working Paper 32, 2009. Available from: http://biostats.bepress.com/upennbiostat/art32 [Accessed on 22 September 2013].
Supporting information Additional supporting information may be found in the online version of this article at the publisher’s web site.
3404 Copyright © 2014 John Wiley & Sons, Ltd.
Statist. Med. 2014, 33 3398–3404
Accepted 10 April 2014
Published online 8 May 2014 in Wiley Online Library
(wileyonlinelibrary.com) DOI: 10.1002/sim.6194
A note on implementation of decaying product correlation structures for quasi-least squares‡ Justine Shultsa*† and Matthew W. Guerrab This note implements an unstructured decaying product matrix via the quasi-least squares approach for estimation of the correlation parameters in the framework of generalized estimating equations. The structure we consider is fairly general without requiring the large number of parameters that are involved in a fully unstructured matrix. It is straightforward to show that the quasi-least squares estimators of the correlation parameters yield feasible values for the unstructured decaying product structure. Furthermore, subject to conditions that are easily checked, the quasi-least squares estimators are valid for longitudinal Bernoulli data. We demonstrate implementation of the structure in a longitudinal clinical trial with both a continuous and binary outcome variable. Copyright © 2014 John Wiley & Sons, Ltd. Keywords:
decaying product; feasible; generalized estimating equations; longitudinal binary data; quasi-least squares
1. Introduction ( )′ We consider longitudinal studies, in which vectors of outcome measurements Y i = Yi1 , … , Yini and ′ ′ associated covariates xij = (xij1 , … , xijp ) are collected at times Ti = (ti1 , … , tini ) on participant i (i = 1, … , m). The expected value and variance of Yij on subject i can be expressed using a generalized linear model: ( ) ( ) ( ) ( ) (1) E Yij = g−1 xij′ 𝛽 = 𝜇ij and Var Yij = 𝜙h 𝜇ij , respectively. For Bernoulli variables, we assume the logistic model, with g−1 (𝛾) = exp(𝛾)∕(1 + exp(𝛾)) = pij , 𝜙 = 1, and h(𝜇ij ) = pij (1−pij ). For continuous variables, we assume the linear model, with g−1 (𝛾) = 𝛾 and h(𝜇ij ) = 1. We also assume that observations on different subjects are independent, while the repeated measurements are correlated within subjects. This note models the correlation using quasi-least squares (QLS) for unbalanced data [1–3] to implement a decaying product correlation structure. QLS is a two-stage computational approach for estimation of the correlation parameters in the framework of generalized estimating equations (GEE). QLS extends GEE by providing a straightforward approach to implement correlation structures such as the decaying product structure, which previously were not available for GEE. For some structures [1], QLS guarantees an estimated correlation matrix that is feasible (positive-definite), while there is no similar guarantee for the corresponding GEE moment-based estimator. As a result, QLS is sometimes suggested as an alternative approach, should the standard implementation of GEE yield an infeasible estimate of the correlation structure that results in a failure to converge [4]. a Department of Biostatistics and Epidemiology, University of Pennsylvania School of Medicine, Philadelphia, b Division of Biometrics III, Office of Biostatistics, OTS, CDER, FDA, Silver Spring, MD, U.S.A.
PA, U.S.A.
*Correspondence
3398
to: Justine Shults, Department of Biostatistics and Epidemiology, University of Pennsylvania School of Medicine, 423 Guardian Drive, 635 Blockley Hall, Philadelphia, PA 19104, U.S.A. † E-mail: [email protected] ‡ The views expressed in this paper are the authors’ personal views and not necessarily those of the US Food and Drug Administration.
Copyright © 2014 John Wiley & Sons, Ltd.
Statist. Med. 2014, 33 3398–3404
J. SHULTS AND M. W. GUERRA
∏ Decaying product correlations of the form Corr(Yij , Yik ) = k−1 w=j Ciw,w+1 (𝛼) (j < k) are plausible to describe the pattern of association among the repeated measurements on a participant in a longitudinal study. Specifying Cij,j+1 (𝛼) = 𝛼 yields the first order autoregressive or AR(1) structure, for which Corr(Yij , Yik ) = 𝛼 k−j . Setting Cij,j+1 (𝛼) = 𝛼 tij+1 −tij yields the Markov structure, for which Corr(Yij , Yik ) = 𝛼 tik −tij . Setting Cij,j+1 (𝛼) = 𝛼j,j+1 yields the unstructured decaying product, or first-order antedependent ∏k−1 AD(1) structure, for which corr(Yij , Yik ) = w=j 𝛼w,w+1 . For the AD(1) structure 𝛼 = (𝛼1,2 , … , 𝛼N−1,N )′ , where N is the number of pre-planned fixed measurement occasions. [1, 2] implemented the AR(1) and Markov structures for QLS analysis of unbalanced data, while [5] implemented the AR(1), Markov, and AD(1) structures for first-order maximum likelihood (MARK1ML) analysis of longitudinal Bernoulli data. As explained in Section 3.2 of [5], the AR(1) structure is often applied for analysis of longitudinal data that are equally spaced in time, while the more general Markov structure is utilized for unequally spaced measurements. The AD(1) structure can be applied when the number of measurements varies between subjects, but might be most plausible when the subjects have a common set of timings with few intermittent missed visits, that is, most subjects who miss a visit do not return to the study [5]. This note implements the AD(1) structure for QLS. Utilizing the tri-diagonal form of the inverse of the AD(1) structure, it is straightforward to show that the QLS estimators of 𝛼j,j+1 are of the same form as those that were previously derived for the AR(1) structure [1, 2], but restricted to measurement occasions j and j + 1 only. As a result, the estimators are feasible, which means that the QLS estimates always yield a positive definite correlation matrix. Furthermore, subject to conditions that are easily checked [6], the estimators will be valid for longitudinal Bernoulli data, which means that there will be a multivariate distribution with the estimated marginal means and correlation structure. This is important because, while a multivariate (normal) distribution can always be constructed for a given mean and correlation structure for continuous random variables, there is no similar guarantee for Bernoulli variables [7–9]. The AD(1) structure is therefore especially useful for the analysis of longitudinal Bernoulli data because its validity is easily checked and it is a flexible structure with fewer correlation parameters than in a fully unstructured matrix, a feature that is also beneficial for continuous outcomes [10].
2. Implementation of the AD(1) structure for quasi-least squares The following algorithm can be used to implement the decaying product structure for QLS. For simplicity, our presentation assumes a common set of measurement times and that subjects have no intermittent missed visits. However, in the Supporting Materials for the online version of this article (Online Appendix B), we generalize the approach to allow some visits to be missed intermittently. Stage one of QLS: First, select starting values for the stage one estimates 𝛼̃ j−1,j . In our experience, starting values of zero are convenient and perform well, for example, 𝛼̃ j−1,j = 0 ( j = 2, …, N; N = ( )′ ̃ QONE = 𝛼̃ 1,2 , … , 𝛼̃ N−1,N . Next, alternate between Steps A and B: max{ni }). Let 𝜶 Step A: Update the estimate of regression parameter 𝜷 by solving the GEE estimating equation [11] for 𝜷: m ∑
−1∕2 −1 Ri
D′i Ai
) ( ) −1∕2 ( 𝛼̃ QONE Ai Y i − Ui (𝜷) = 0,
(2)
i=1
( )′ ∏ where Ui (𝜷) = 𝜇i1 , 𝜇i2 , · · · , 𝜇ini , Di = 𝜕Ui (𝜷)i ∕𝜕𝜷, Ri [j, k] = k−1 w=j 𝛼w,w+1 (j < k), Ri [k, j] = Ri [j, k], )) ( ( and Ai = diag h(𝜇i1 ), · · · , h 𝜇inii . Step B: Update the estimate of the correlation by obtaining a solution (for 𝜶) to the stage one estimating ̂ for 𝜷: equation evaluated at the current estimate 𝜷 ( ) 𝜕R−1 (𝜶) ( ) i ̂ ̂ = 0 (j = 2, … , N), Zi′ 𝜷 Zi 𝜷 𝜕𝛼 j−1,j i=1
m ∑
where
(
) ( )′ Y i − Ui (𝜷) = zi1 , … , zini
3399
−1∕2
Zi (𝜷) = Ai
(3)
is the vector of Pearson residuals on subject i. Copyright © 2014 John Wiley & Sons, Ltd.
Statist. Med. 2014, 33 3398–3404
J. SHULTS AND M. W. GUERRA
Stage two of QLS: After convergence in the parameters in stage one, we next obtain a consistent estimate of 𝜶 by solving the stage two estimating equation for 𝜶: m ∑
( trace
i=1
𝜕R−1 (𝜶) || i Ri (𝜶) | 𝜕𝛼j−1,j ||𝜶=𝛼̃ QONE
) = 0 (j = 2, … , N).
(4)
( )′ ̂1,2 , … , 𝛼 ̂N−1,N . We then solve the GEE estimating We refer to the solution of (4) as 𝛼 ̂QLS = 𝛼 ( ) equation (2) at R−1 𝛼 ̂QLS , in order to obtain the final estimate of 𝜷. i To simplify (3) and (4), we utilize the simple tri-diagonal form for the inverse of Ri (𝜶) of ) the inverse of for an ni × ni (ni ∈ {2, … , N}) decaying product structure. The elements ( −1 2 [j, j] = ni × ni decaying product AD(1) matrix Ri (𝜶) are [12]: Ri [1, 1] = 1∕ 1 − 𝛼1,2 ; R−1 i ( ) (( )( )) 2 2 2 2 −1 1 − 𝛼j−1,j 𝛼j,j+1 ∕ 1 − 𝛼j−1,j 1 − 𝛼j,j+1 (j = 2, … , ni − 1 for ni > 2); Ri [j, j + 1] = ) ) ( ( 2 2 −𝛼j,j+1 ∕ 1 − 𝛼j,j+1 (j = 1, … , ni − 1); R−1 ; and R−1 [n , n ] = 1∕ 1 − 𝛼 [j, j′ ] = 0 for |j − j′ | > 1. i i i i ni −1,ni The derivative of the inverse of ni × ni (ni ⩾ 2) decaying product matrix Ri (𝜶) with respect to 𝛼j−1,j is then easily shown to equal an ni × ni matrix with the following 2 × 2 sub-matrix in rows j − 1 to j and columns j − 1 to j: 𝜕R−1 (𝜶) [ i 𝜕𝛼j−1,j
) ( 2 ⎛ ⎞ − 1 + 𝛼 2𝛼 j−1,j 1 j−1,j ⎟ ⎜ ( ) j − 1 ∶ j, j − 1 ∶ j = ( )2 ⎜ 2 ⎟ 2𝛼j−1,j 2 ⎝ − 1 + 𝛼j−1,j ⎠ 1 − 𝛼j−1,j ]
(j = 2, … , ni )
(5)
and, for ni > 2, all other elements equal to 0. The solutions to the stage one and stage two estimating equations (3) and (4) can then be directly obtained as follows: √ ) m ( m ( m ( )2 ∑ )2 ∑ ∑ ̂z2ij + ̂z2ij−1 − ̂zij + ̂zij−1 ̂zij − ̂zij−1 𝛼 ̃j−1,j =
i=1
i=1
∑ m
2
i=1
( j = 2, … , N)
(6)
̂ziĵzij−1
i=1
and 𝛼 ̂j−1,j =
2𝛼̃ j−1,j 2 1 + 𝛼̃ j−1,j
( j = 2, … , N).
(7)
The summations in (6) are taken over subjects who have measurements at both visits j and j - 1, so that 𝛼 ̂1,2 (for example) is based on ̂zi1 and ̂zi2 from subjects who were present at visits 1 and 2. The estimators given in (6) and (7) are of the same form as the stage one and two estimators for the AR(1) structure [1, 2], respectively, for a subject with measurements at visits j − 1 and j only. Therefore, the proofs that were provided earlier [1, 2] are applicable to show that both −1 < 𝛼 ̃j−1,j < 1 and −1 < 𝛼 ̂j−1,j < 1 ∀j = 2, … , N; as a result, the estimated AD(1) matrix is positive definite, so that the QLS estimators are feasible for this structure. In addition, a valid multivariate distribution can be constructed ̂ and 𝜶 ̂ QLS as long as the estimates satisfy the Prentice using model (6) in [5] for the given QLS estimates 𝜷 constraints [6] for the consecutive correlations: ̂j−1,j ⩽ Upperi ( j − 1, j), Loweri ( j − 1, j) ⩽ 𝛼
(i = 1, … , m ; j = 2, … , ni )
(8)
3400
{ ( )1∕2 ( )−1∕2 } where Loweri (j − 1, j) = max − w , Upperi (j − 1, j) = min ̂ ij−1 w ̂ ij ,− w ̂ ij−1 w ̂ ij {( )1∕2 ( )−1∕2 } ( )−1 wij , w ̂ ij−1 ∕̂ wij pij 1 − ̂ w ̂ ij−1 ∕̂ pij . In general, the QLS estimators for all , and w ̂ ij = ̂ decaying product structures (Markov, AR(1), AD(1)) will be valid for Bernoulli data as long as the estimated consecutive correlations satisfy the standard Prentice constraints (8) for a bivariate Bernoulli distribution [5]. For other structures, for example, the exchangeable structure with Corr(Yij , Yik ) = 𝛼 for all i, j, k, satisfaction of the constraints (8) must be checked for all pairs. Copyright © 2014 John Wiley & Sons, Ltd.
Statist. Med. 2014, 33 3398–3404
J. SHULTS AND M. W. GUERRA
2.1. Demonstration We next implement the AD(1) structure for a QLS analysis of a longitudinal trial to compare Effexor with lithium in the treatment of bipolar type II depression. Approximately 30% of measurements were randomly dropped in an earlier demonstration, but here we consider the full data set from the Effexor study that was described previously [3, Section 9.3]. We implement model (1) with logit and identity links for a binary and continuous outcome, respectively. In addition, xij′ 𝛽 = 𝛽0 + 𝛽1 I(Effexor) + 𝛽2 visit + 𝛽3 I(Effexor) × visit
(9)
where I(Effexor) takes value 1 or 0 for the 26 patients treated with Effexor or lithium, respectively, and visit = 1, … , ni (ni ∈ {2, … , 8}). There were no subjects with intermittent missing measurements and the mean (median) number of measurements was 6.4 (8.0). If the regression parameter 𝛽3 that is associated with the Effexor by visit interaction differs significantly from zero, this will indicate that the mean change over time in HAMD-17 depression scores (continuous outcome) or log odds of depression (Bernoulli outcome) differs significantly between the two treatment groups. Depression was defined as having a HAMD-17 score > 7. Table I displays the results of the analysis for the HAMD-17 scores, while Table II shows the results for the binary outcome. Several working correlation structures were implemented in the analysis. Results for the fully unstructured matrix are not shown because convergence was not achieved for this structure, for either the binary or the continuous outcome. The DBAR, CONE, and CTWO Rotnizky-Jewell criteria [13] that were assessed in [14] were used to aid in the selection of a working correlation structure. Values of (DBAR, CONE, CTWO) that are closer to (0, 1, 1) indicate better fit. For the continuous outcome, the DBAR criterion selected the AD(1) structure, while the CONE and CTWO criteria selected the AR(1) structure. However, the fit was similar and good for both the AD(1) and AR(1) structures, worse for the exchangeable structure, and very poor for the identity structure. For
Table I. Estimated values of the regression and correlation parameters in the analysis of data from a clinical trial to compare Effexor versus lithium with respect to mean change in HAMD-17 depression scores. Structure AR(1) EXC AD(1) IND
𝛽̂0
𝛽̂1
𝛽̂2
𝛽̂3
𝛼 ̂QLS
pval
DBAR
CONE
CTWO
22.57 21.99 21.93 21.71
−0.33 −1.04 −0.53 −0.65
−1.21 −1.30 −1.06 −1.17
−0.69 −0.54 −0.67 −0.67
0.76 0.62 ** 0
0.05 0.11 0.07 0.11
0.31 0.95 0.29 10.19
0.95 1.50 0.99 3.14
1.21 2.95 1.28 15.46
Results are shown for the AR(1), exchangeable (EXC), AD(1), and identity structures. The pvalue (pval) for the two-sided test of Ho ∶ 𝛽3 = 0 and the CONE, CTWO, and DBAR criteria for selection of a working correlation are also provided. ̂1,2 = 0.68, 𝛼 ̂2,3 = 0.81, 𝛼 ̂3,4 = 0.62, 𝛼 ̂4,5 = 0.69, 𝛼 ̂5,6 = 0.82, 𝛼 ̂6,7 = 0.80, 𝛼 ̂7,8 = 0.90 ** 𝛼
Table II. Estimated values of the regression and correlation parameters in the analysis of data from a clinical trial to compare Effexor versus lithium with respect to change in the log odds of depression. Structure AR(1) EXC AD(1) IND
𝛽̂0
𝛽̂1
𝛽̂2
𝛽̂3
𝛼 ̂QLS
Lowers
Upper
pval
DBAR
CONE
CTWO
0.47 0.50 0.48 0.43
0.42 −0.28 0.55 0.14
−0.22 −0.28 −0.23 −0.22
−0.33 −0.31 −0.36 −0.27
0.55 0.41 * 0
−0.04 −0.02 ** ***
0.76 0.13 ** ***
0.10 0.30 0.07 0.18
0.42 5.69 0.53 5.02
1.30 2.42 1.37 2.74
2.02 9.52 2.27 9.51
Copyright © 2014 John Wiley & Sons, Ltd.
Statist. Med. 2014, 33 3398–3404
3401
Results are shown for the AR(1), exchangeable (EXC), AD(1), and identity structures. The p-value (pval) for the two-sided test of Ho ∶ 𝛽3 = 0 and the CONE, CTWO, and DBAR criteria for selection of a working correlation are ̂QLS are also provided. also provided. The lower (Lower) and upper (Upper) Prentice constraints (8) for 𝛼 ̂1,2 = 0.60(−0.64, 0.74); 𝛼 ̂2,3 = 0.68(−0.36, 0.74); 𝛼 ̂3,4 = 0.54 (-0.2,0.74); 𝛼 ̂4,5 = 0.54(−0.20, 0.74); 𝛼 ̂5,6 = *𝛼 0.29(−0.11, 0.74); 𝛼 ̂6,7 = 0.51(−0.06, 0.74); 𝛼 ̂7,8 = 0.59(−0.03, 0.74) ** The lower and upper Prentice constraints (8) are provided in parentheses following each aforementioned estimate. *** The Prentice constraints (8) are always satisfied for the identity structure.
J. SHULTS AND M. W. GUERRA
example, DBAR was 0.29 and 0.31 for the AD(1) and AR(1) structures, versus 0.95 and 10.19 for the exchangeable and identity structures, respectively. Assessment of the sample variogram [15, p. 49] also suggested that there is a decay in the intra-subject correlations that is not compatible with an identity or an exchangeable working structure. The sample variogram is described in the Supporting Materials for the online version of this article (Online Appendix A). For the binary outcome, all criteria selected the AR(1) structure. However, the fit was similar for the AR(1) and AD(1) structures and was much worse for both the exchangeable and identity structures. For example, DBAR was 0.42 and 0.53 for the AR(1) and AD(1) structures versus 5.69 and 5.02 for the exchangeable and identity structures, respectively. It is straightforward to show [9, Table 1] that the upper Prentice constraint for model (9) for both the AR(1) and AD(1) structures is { ( ) ( )} min exp −|𝛽̂2 |∕2 , exp −|𝛽̂2 + 𝛽̂3 |∕2 , while the lower Prentice constraint is always negative. The Prentice constraints were not violated for the AR(1) and AD(1) structures, but the estimate 𝛼 ̂QLS = 0.41 clearly violated the estimated Prentice constraints (−0.02, 0.13) for the exchangeable structure. The exchangeable structure should therefore be ruled out according to the rule out criterion [14]. The Prentice constraints will never be violated for the identity structure because the lower bound for the constraint is always negative, while the upper bound is positive. However, the poor fit of the identity structure suggests that it should also be ruled out in favor of the AR(1) or AD(1) structures, which had better (and similar) fit. The results of the analysis were similar for all structures, although the results for the AR(1) and AD(1) structures were most suggestive of a beneficial impact of the intervention. For example, if the criterion for statistical significance was a p-value ⩽ 0.1, then the results for both the continuous and binary outcomes would be significant only for the AD(1) and AR(1) structures. However, it is important to acknowledge that we implemented sandwich-type standard errors that can be underestimated in small samples [16]. Work is needed to assess small sample corrections to sandwich-type estimators of covariance for QLS, similar to assessments made for GEE [17]. The results in Tables 1 and 2 were based on a full data set with no intermittent missed measurements. However, in the Supporting Materials for the online version of this article (Online Appendix B), we present the results for the same data set, but with some randomly dropped measurements and intermittently missed visits. The Supporting Materials for the online version of this article (Online Appendix C) present the results of a small simulation study to assess the implementation of the AD(1) structure for QLS when model (1) is correctly specified and several working structures (Markov, exchangeable, unstructured, and AD(1)) are implemented in the analysis. We considered settings for which the Markov structure was the only true structure, and for which the AD(1) and unstructured matrices were the only true structures. We also assessed the conditions under which violation of the Prentice constraints (8) is more likely to occur. The simulations demonstrated that correctly specifying the working correlation structure, whether Markov or AD(1), resulted in greater small sample efficiency in estimation of the regression parameters. In addition, failure to converge was not an issue when fitting models to the longitudinal Bernoulli data that we simulated using the algorithms in [18, 19]. The simulations also suggested that violation of the Prentice constraints is not likely for a correctly specified working structure, unless the true correlation is very close to a boundary value for the constraints; even then, the violation is not likely to be severe (as assessed by the distance between the estimated correlation parameter and the estimated closest boundary value), especially for increasingly larger sample sizes. However, if the correlation structure is incorrectly specified, then there can be a severe violation of constraints that is not alleviated by increasing sample size. For example, we considered settings for which 𝛼̂ for the incorrectly specified exchangeable structure converges in probability to a value that exceeds the upper boundary value of the Prentice constraints; for these scenarios, the percentage of simulation runs that resulted in a violation of constraints for the incorrectly specified exchangeable structure increased to 100% when the group size increased to m = 60.
3. Conclusion
3402
This note implemented the AD(1) structure for QLS and demonstrated the results in analysis of continuous depression scores and Bernoulli depression status in a randomized clinical trial. Expanding the list of correlation structures available for GEE and QLS allows for more careful modeling of the correlation structure for GEE, which is helpful in avoiding the serious loss in efficiency that can result from Copyright © 2014 John Wiley & Sons, Ltd.
Statist. Med. 2014, 33 3398–3404
J. SHULTS AND M. W. GUERRA
incorrect application of an identity structure [20–23]. In addition, unlike the estimators of the regression parameters, the GEE and QLS estimators of the association parameters may fail to be consistent when the true correlation structure is misspecified; in this situation, Sutradhar and Das [21] suggest that the sandwich-based covariance matrix should be evaluated at the limiting values of the association parameters. As a result, application of the sandwich covariance matrix to correct the standard errors for potential misspecification of the true structure can also require careful assessment of the correlation structure. The AD(1) structure is more flexible than the AR(1) or exchangeable structures that assume a constant value for the consecutive and all pairwise correlations, respectively. The Markov structure that depends on the spacing of measurements might be more appropriate when the temporal spacing of measurements is highly irregular, either due to study design (as in [24]) or to many intermittently missing visits. The assumption of a decaying product pattern in the correlations is not as flexible as the fully unstructured matrix, which can be implemented for balanced data (all ni = n) for QLS [3, Section 5.6]. However, the increased flexibility of the fully unstructured matrix can come at a cost of problems with respect to convergence, which was the case for both outcomes of this study, when implementation of the unstructured matrix was attempted with GEE. The analysis was conducted using the xtqls command [4] in Stata 13.0 that we updated to fit the AD(1) structure. We also updated the qls function utilized in the qlspack package [25] to allow for implementation of the AD(1) structure in R. Stata code to recreate the analysis results in Section 2.1 and Online Appendices A and B is provided in the Supporting Materials for the online version of the article (Online Appendix D). R code to recreate the analysis results in Section 2.1 is provided in Online Appendix E.
Acknowledgements The authors are very grateful to the associate editor and two anonymous reviewers for their careful review and thoughtful suggestions, which allowed us to make considerable improvements to the manuscript.
References
Copyright © 2014 John Wiley & Sons, Ltd.
Statist. Med. 2014, 33 3398–3404
3403
1. Shults J, Chaganty NR. Analysis of serially correlated data using quasi-least squares. Biometrics 1998; 54(4): 1622–1630. 2. Chaganty NR, Shults J. On eliminating the asymptotic bias in the quasi-least squares estimate of the correlation parameter. Journal of Statistical Planning and Inference 1999; 76:145–161. 3. Shults J, Hilbe J. Quasi-Least Squares Regression. Chapman & Hall/CRC Press Monographs on Statistics & Applied Probability: Boca Raton, 2014. 4. Shults J, Ratcliffe S, Leonard M. Improved generalized estimating equation analysis via xtqls for implementation of quasileast squares in Stata. Stata Journal 2007; 7:147–166. 5. Guerra MW, Shults J, Amsterdam J, Ten-Have T. The analysis of binary longitudinal data with time-dependent covariates. Statistics in Medicine 2012; 10:931–948. 6. Prentice RL. Correlated binary regression with covariates specific to each binary observation. Biometrics 1988; 44: 1033–1048. 7. Ziegler A, Vens M. Generalized estimating equations: notes on the choice of the working correlation matrix. Methods of Information in Medicine 2010; 5:421–425. 8. Molenberghs G. Editorial. “Generalized estimating equations: notes on the choice of the working correlation matrix”. Methods of Information in Medicine 2010; 49:419–420. 9. Shults J. Discussion of “Generalized estimating equations: notes on the choice of the working correlation matrix” continued. Letter to the editor. Methods of Information in Medicine 2011; 50:96–99. 10. Westgate PM. A bias correction for covariance estimators to improve inference with generalized estimating equations that use an unstructured correlation matrix. Statistics in Medicine 2012; 35(13):2850–2858. 11. Liang KY, Zeger SL. Longitudinal data analysis using generalized linear models. Biometrika 1986; 73:13–22. 12. Zimmerman DL, Nu˜nez-Antón VA. Antedependence Models for Longitudinal Data. Chapman and Hall/CRC Press: Boca Raton, 2012. 13. Rotnitzky A, Jewell NP. Hypothesis testing of regression parameters in semiparametric generalized linear models for cluster correlated data. Biometrika 1990; 77:485–497. 14. Shults J, Sun W, Tu X, Kim H, Amsterdam J, Hilbe JM, Ten-Have T. A comparison of several approaches for choosing between working correlation structures in generalized estimating equation analysis of longitudinal binary data. Statistics in Medicine 2009; 28:2338–2355. 15. Diggle P, Heagerty P, Liang KY, Zeger S. Analysis of Longitudinal Data. Oxford University Press: New York, 2002. 16. Mancl LA, DeRouen TA. A covariance estimator for GEE with improved small-sample properties. Biometrics 2001; 57(1):126–134. 17. Teerenstra S, Bing L, Preisser JS, Achterberg T, Borm G. Sample size considerations for GEE analyses of three-level cluster randomized trials. Biometrics 2010; 66(4):1230–1237.
J. SHULTS AND M. W. GUERRA 18. Guerra MW, Shults J. A note on the simulation of overdispersed random variables with specified marginal means and product correlations. The American Statistician. In Press; Available from: http://www.tandfonline.com/doi/ [Accessed on 18 Febuary 2014]. 19. Preisser JS, Qaqish BF. A comparison of methods for simulating correlated binary variables with specified marginal means and correlations. Journal of Statistical Computation and Simulation. In Press; Available from: http://www.tandfonline. com/ [Accessed on 18 February 2014]. DOI: abs/10.1080/00949655.2013.818148. 20. Fitzmaurice G. A caveat concerning independence estimating equations with multivariate binary data. Biometrics 1995; 51(1):309–317. 21. Sutradhar B, Das K. On the efficiency of regression estimators in generalised linear models for longitudinal data. Biometrika 1999; 86(2):459–465. 22. Sutradhar B, Das K. On the accuracy of efficiency of estimating equation approach. Biometrics 2000; 56(2):622–625. 23. Wang Y, Carey V. Working correlation structure misspecification, estimation and covariate design: implications for generalised estimating equations performance. Biometrika 2003; 90(1):29–41. 24. Preisser J, Arcury T, Quandt S. Detecting patterns of occupational illness clustering with alternating logistic regressions applied to longitudinal data. American Journal of Epidemiology 2003; 158(5):495–501. 25. Xie J, Shults J. Implementation of quasi-least squares With the R package qlspack. UPenn Biostatistics Working Papers. Working Paper 32, 2009. Available from: http://biostats.bepress.com/upennbiostat/art32 [Accessed on 22 September 2013].
Supporting information Additional supporting information may be found in the online version of this article at the publisher’s web site.
3404 Copyright © 2014 John Wiley & Sons, Ltd.
Statist. Med. 2014, 33 3398–3404
