Authors: Lin GE, Justin X. TU, Hui ZHANG, Hongyue WANG, Hua HE, Douglas GUNZLER
Categories: Biostatistics in Psychiatry (35), R, SAS, binary variables, correlated outcomes, generalized linear mixed-effects models, latent variable models, weighted generalized estimating equations, 二分类变量, 加权广义估计方程, 广义线性混合效应 模型, 潜变量模型, 相关结果
Source: Shanghai Archives of Psychiatry
Authors: Lin GE, Justin X. TU, Hui ZHANG, Hongyue WANG, Hua HE, Douglas GUNZLER
Longitudinal studies are used in mental health research and services studies. The dominant approaches for longitudinal data analysis are the generalized linear mixed-effects models (GLMM) and the weighted generalized estimating equations (WGEE). Although both classes of models have been extensively published and widely applied, differences between and limitations about these methods are not clearly delineated and well documented. Unfortunately, some of the differences and limitations carry significant implications for reporting, comparing and interpreting research findings.
In this report, we review both major approaches for longitudinal data analysis and highlight their similarities and major differences. We focus on comparison of the two classes of models in terms of model assumptions, model parameter interpretation, applicability and limitations, using both real and simulated data. We discuss caveats and cautions when applying the two different approaches to real study data.
Longitudinal study designs have become increasingly popular in research and practice across all disciplines. Such designs capture both between-individual differences and within-subject dynamics, providing opportunities to study complex biological, psychological and behavioral changes over time such as causal treatment effects and mechanisms of change ^[1,2]^. Since longitudinal study designs create serial correlations over repeated assessments from same subjects, traditional statistical methods for cross-sectional data analysis such as linear and logistic regression do not apply. In addition, since longitudinal studies are typically of long duration, missing data is common. Specialized models and methods must be used to address the two major issues.
The two dominant approaches for longitudinal data analysis are the generalized linear mixed-effects model (GLMM) and weighted generalized estimating equations (WGEE)^[1]^. Both methods are derived from the same class of models for cross-sectional data, the generalized linear models (GLM). Because different techniques are used to extend the GLM to longitudinal data, the GLMM and WGEE are quite different and some of the differences carry significant implications for their applicability and interpretation of study findings. Despite the existence of an extremely large body of literature discussing the development and application of the two approaches, practitioners still often face many difficult questions when choosing and applying such models to real study data.
For example, which approach is applied given data from a study? Do the two approaches yield identical estimates and/or inferences? If not, how should one approach and interpret such differences? What are the pros and cons associated with each approach? Although some questions have well-documented answers in the literatures, others have only been approached recently and still await answers.
In this report, we first give an overview of the approaches and then discuss major differences between the two classes of models. Unlike the literature on the discussion of the two methods, we focus on their practical implications, which we think provide useful guidance for practitioners for selecting right approaches for their studies and effectively addressing their study questions.
Since both the GLMM and WGEE are extensions of the GLM, we start with a brief overview of the latter.
Consider a sample of n subjects and let Yi (Xi) denote a continuous response (a vector of explanatory variables). The classic linear model is given
(1)
where N(μ, σ^2^) denotes a normal distribution with mean μ and σ^2^ variance. The linear model is widely used in research and practice. One major limitation is that it only applies to continuous response Yi. The generalized linear models (GLM) extend the classic linear model to non-continuous response such as binary.
To express the GLM, we first rewrite the linear regression in (1) as
(2)
where yi | Xi denotes the conditional distribution of Yi given Xi and E(yi |xi) denotes the conditional mean of Yi given Xi. By replacing the normal in (2) with other distributions appropriate for the type of response, we obtain the class of GLM:
(3)
where f(μ) denotes some distribution with mean μ and g(μ) is a function of μ. Since g(μ) links the mean to the explanatory variables, g(μ) is called the link function.
The specification of f(μ) and g(μ) depend on the type of response Yi. For a binary Yi, f(μ) is the Bernoulli distribution and g(μ) is often set as the function logit .
The resulting GLM is the logistic regression. For a count Yi, a popular choice for f(μ) is the Poisson distribution and g(μ) is the log function, g(μ) =log(μ). Another choice for f(μ) is the negative binomial. A major limitation of Poisson is that its variance is the same as its mean. Count responses arising in many real studies often have variances larger than means, a phenomenon known as “overdispersion” ^[1]^. The negative binomial is similar to the Poisson, but unlike the Poisson, allows for overdispersed count responses ^[1]^.
Inference for GLM can be based on maximum likelihood (ML) or estimating equations (EE). The classic ML provides most efficient estimates, if the response Yi follows the specified distribution such as the normal in the linear regression in (1). In many studies, it may be difficult to specify the right distribution, in which case the ML will yield biased estimates if the specified distribution does not match the data distribution. The modern alternative EE uses an approach for inference that does not require specification of a mathematical distribution for Yi, thereby providing valid inference for a wider class of data distribution. Since no distribution is required under EE, we may also express the GLM in this case
(4)
or simply
(5)
or equivalently
(6)
where h = g^–1^ is the inverse of g(μ). When specified without the distribution component, (4), (5) or (6) are also called the semi-parametric GLM. In comparison, (3) is called the parametric GLM.
Example 1. In a suicide study, we may model number of suicide attempts, Yi, a count response, using a parametric GLM:
(7)
where Xi denotes a set of explanatory variables such as age, physical problems and history of depression. If concerned about overdispersion, such as indicated empirically by a much larger variance as compared to the mean, the semi-parametric GLM below may be used
Unlike the parametric model in (7), the semi-parametric GLM above provides valid inference regardless of presence of over-dispersion.
Consider a longitudinal study with T time points and let Yit and Xit denote the same response and predictors/covariates as in the cross-sectional setting, but with t indicating their dependence on the time of assessment (1<i<n,1<t<T). For each time t, we can apply the GLM in (3) to model the regression relationship between Yit and Xit at each point
We can then get estimates of βt for each time point t. However, it is difficult to interpret different βt across the different time points. Moreover, it is technically challenging to combine estimates of βt to test hypotheses concerning temporal trends because of interdependence between such estimates.
The WGEE addresses the aforementioned difficulties by using a single estimate β to model changes over time based on multiple assessment times^[1,3]^. Since the WGEE estimates β using a set of equations that do not reply on assumed distribution f(μit) in (8), the first part of the GLM can be removed and the resulting model becomes
(8)
By comparing the WGEE above with the model in (4), it is seen that the WGEE is an extension of the semi-parametric GLM to longitudinal data. The key difference between (4) and (8) is that (8) is not simply an application of GLM to each of the time points, but rather an extension of the model in (4) to provide a single parameter vector β for easy interpretation and estimate this parameter vector by using data from all time points and accounting for correlations between the repeated assessments. Like the semi-parametric GLM, the WGEE provides valid inference for a wider class of data distributions.
The GLMM extends the GLM to longitudinal data analysis using a completely different approach.
Consider again a longitudinal study with T time points and let Yit and Xit denote the same response and predictors/covariates as in the WGEE above. The GLMM is specified
(9)
where N(μ, N (μ, Σ) denotes a multivariate normal with mean μ and variance Σ, Zit is a vector of predictors/covariates (often set equal to Xit), and g(μ) is the appropriate link function for the type of response Yit. The vector of latent variables, bi, is called the random effects, denoting individual differences from the population mean bi, which is known as the fixed effects. Although β̃ is typically assumed to follow a multivariate normal as in (3), other types of distributions may also be considered.
Unlike the WGEE, the GLMM accommodates correlated responses Yit by directly modeling their joint distribution. Since multivariate distributions are extremely complex except for the multivariate normal, latent variables bi are generally employed to model the correlated responses. Thus, although Yit is still modeled for each time point t, by including the random effect b, in the specification of the conditional distribution of Yit given bi (Xit and Zit), the GLMM in (9) allows the resulting Yit ‘s to be correlated (conditional on Xit and Zit only). This approach allows one to specify multivariate distributions using familiar univariate distributions such as the Bernoulli (for binary responses) and the Poisson (for count responses).
Although both WGEE and GLMM are extensions of the GLM, the two approaches are quite different because of the way the extensions are accomplished. In this section, we discuss such key differences.
A fundamental difference between the WGEE and GLMM is in interpretation of model parameters. As noted earlier, the WGEE is an extension of the semi-parametric GLM, while the GLMM is an extension of the parametric GLM. Although the parameter vector β has the same interpretation between the parametric (1) and semi-parametric (4) GLM, the parameter vectors β in the WGEE (8) and β in the GLMM (9) are generally different except for the linear regression.
Example 2. Consider the model in Example 1, but now assume a longitudinal study with T assessment times. Let Yit and Xit denote the longitudinal versions of Yi and Xi. Again, if Yit at each point is modeled to follow a Poisson, the GLMM has the
(10)
On the other hand, the WGEE for a count response has the form
(11)
The two models look quite the same, except for the additional random effect in (10). So, to compare the two models, we integrate out bi in (10) and calculate the conditional expectation E (yit | xit) to obtain (see Zhang et al., 2012 for derivationμ]
(12)
This conditional expectation is different from the expression in (11) for the WGEE unless
In some special cases, 3 and [3 may be identical except for the intercept term. For example, if z is independent of xit, then
Thus (12) reduces to
(13)
Except for the intercept term, β and i have the same interpretation. The next Example illustrates this special case with data from a real study.
Example 3. The COMBINE study was a multi-site randomized clinical trial with longitudinal follow-up to compare intervention effects between nine pharmacological and/or psychosocial treatments. For illustration purposes, we combined all the 9 treatment conditions and focused on the two drinking outcomes, days of any drinking and days of heavy drinking over the past month, at three visits during treatment at weeks 8, 16 and 26.
Denote the three visits by 1<t<3. To see whether the mean days of any drinking (heavy drinking) varied in the three time periods, we created two dummy variables indicating weeks 16 and 26 as the fixed effect Xi and a random intercept to account for correlations among the repeated assessments for the GLMM
Since Zi= 1 and Xi1t=1 are both constants, they are independent. We fit a WGEE with the same predictors Xit=(Xi1t, Xi2t, Xi3t)^T^ as in the fixed effect of GLMM above. As indicated in Example 2, β̃0 differs from its WGEE counterpart βo by a constant , whereas β̃1 and β̃2 retain the same scales as the corresponding β1 and β2 in the WGEE model.
Shown in Table 1 are the estimates of β and β̃ and associated standard errors, and test statistics and associated p-values for both models. As expected, estimates of β̃1 and β̃2 were quite close to the respective WGEE estimates β1i and β2, but those of ^were much^ smaller (in magnitude) than their WGEE counterparts (0.694 vs. 1.908 for Days of Any Drinking and -0.30 vs. 1.307 for Days of Heavy Drinking).
Since β̃1 and β̃2 (β1 and β2) represent changes relative to the reference level β̃0(β0) from Visit 1 to Visit 2 and from Visit 1 to Visit 3 under GLMM (WGEE), the difference between the estimated β̃0 and β0 implies not only different means at visit 1, but also at visits 2 and 3 for each drinking outcome. For example, for the outcome of Days of Any Drinking, the WGEE estimates indicate a mean of 6. 7, 7.7 and 8.9 days of any drink over the three visits, much lower than 1.9, 2. 3 and 2.7 days of any drink at the corresponding visits estimated by the GLMM. The WGEE estimates are actually identical to sample means of this drinking outcome at each visit, but the GLMM estimates are not, making the latter difficult to interpret.
Example 4. In Example 2, now suppose that yMit is binary and is modeled by a GLMM for binary response as follows
(14)
The corresponding WGEE has the form
(15)
As in Example 2, the two models look the same, except for the additional random effect in (14). To compare the two, we again need to integrate out bi in (14). In general, it is difficult to obtain a closed-form expression for the resulting integral. But, in the special case when Zit is a subset of Xit and Zit =Zit, the integral can be expressed as (see Zhang, Xia et al., 2011 for derivation)^[5]^
(16)
Thus β̃0, is a rescaled β. Even in this special case, it may not easy to interpret β̃, as demonstrated by a real study example next.
Example 5. In a recent study on smoking cessation, 276 subjects participated in multi-component program adapted to seriously mentally ill patients within an outpatient mental health clinic. Out of these subjects, 99 also participated in a formal evaluation, in which interviews were conducted at the point of enrollment (baseline) and again at 3, 6 and 12 months. A primary outcome of the study is the 7-day point prevalence abstinence (defined as no smoking at all in the previous 7 days). We modeled changes of this longitudinal binary abstinence outcome using data from the 99 subjects.
We applied both longitudinal approaches and used three dummy variables to model the rate of 7-day point prevalence abstinence at each of the follow-up times for the WGEE and also included a random intercept for the GLMM model. Thus, the GLMM has the form
while the WGEE is given
Shown in Table 2 are the estimates of β̃ and i and associated standard errors, and test statistics and associated p-values for both models. The estimates were quite different between the two models, although the ratio , a constant for all 1<k<3. This is not a coincidence, but expected, because of (16) and the fact that Ziti in the GLMM model is a subset of Xit.
Also, since H0 : βt =0 under one model implies that H0 holds true for the corresponding coefficient under the other model, same conclusions were obtained regarding statistical significance of the parameters, although the p-values under the two models were slightly different.
The scale difference in the parameters between the two models did have serious implications for the interpretation of estimates from the GLMM. Shown in the last two columns of Table 2 are the rates of 7-day point prevalence abstinence over time estimated by each of the two models. Since the WGEE becomes the logistic regression at each assessment time, it yields estimates identical to the observed rates of 7-day point prevalence abstinence. Estimates from the GLMM were quite different, underestimating the observed rates by over 50%. Even in this simplest case where estimates of GLMM are a scale shift of their GEE counterparts, GLMM estimates are difficult to interpret.
Inference for the WGEE is based on a set of equations. Although it is generally not possible to express solutions (estimates) in closed-form, the equations are readily solved numerically ^[1]^. Inference for the GLMM, however, is much more challenging. Since the likelihood function arising from the GLMM is generally quite complex, involving multidimensional integrals, it is difficult to directly maximize this function except for the special linear models case for continuous outcomes. Different approaches have been proposed to address the computational issues when modeling non-continuous responses using the GLMM.
The approach implemented in most major software such as SAS and R is integral approximation. This approach first approximates the log-likelihood function and then maximizes the approximated function using the Newton and Gauss-Hermite quadratures. However, studies have shown that the approach does not work well as one would expect ^[6]^, especially for modeling binary responses. Below, we highlight some key findings from these studies.
Example 6. Zhang, Lu and Feng et al. (2011) ^[6]^ examined the computational issues by fitting GLMM for simulated binary responses using both SAS and R. They simulated an explanatory variable Xi and the response Yit from the GLMM:
(17)
where β0=β1=1 and τ=0.001 and 2. For τ=0.001, the within-subject correlation was very small and thus negligible, making theYit’s almost independent. For τ=2, the within-subject correlations was about 0.5.
They simulated data from the GLMM in (17) with a sample size n=500, fit the same model to the simulated data using R and SAS, and repeated the process 1,000 times. Shown in Table 3 are the estimates of parameter vector β (under “β0=1” and “β1=1” based on averaging over 1,000 estimates from fitting the model to each simulated data) and associated standard errors (under S.E. × 10 based on sample standard deviations from the Monte Carlo 1,000 estimates of β) from fitting SAS NLMIXED procedure and R lme4 function. For t=0.001, the estimates were quite similar between the two, but for t=2, the SAS NLMIXED procedure provided more accurate estimates. There were more pronounced differences in the standard errors between the two, with the R lme4 consistently yielding lower standard errors than the SAS NILMIXED. Such differences play a significant role in hypothesis testing.
Shown in Table 3 under “Type I error” are the type I error rates for testing the hypothesis, H0 : β1 = 1 vs. Ha : β1 ≠ 1 based on the Wald statistic from the SAS and R procedures. The Wald statistic is where is the standard error for the estimate. The type I error rate was calculated
where q1,0.95 is the 95th percentile of χ^2^1 and T^(m)^W denotes the Wald statistic for testing the hypothesis based on m th simulated data (1 ≤ m ≤1,000). Since the null hypothesis H0 : β1 = 1 is true, the type I error rate α̂ should be close to the nominal value a=0.05, if the SAS and R procedures provided correct inference. Although the SAS NLMIXED did yield type I error rates close to the nominal value, R lme4’s estimates were inflated, especially for the case with higher within- subject correlation t=2.
In addition to the SAS and R procedures, Zhang, Lu and Feng et al. (2011) ^[6]^ also considered other procedures in SAS and R and found that none provided correct estimates. Their conclusion was that the SAS NLMIXED procedure provided more accurate estimates and type I error rates. However, more recent studies by Chen, Knox, Arora et al. (2016) ^[7]^ and Chen, Lu, Arora et al. (2016)^[8]^ show that this SAS procedure did not provide correct estimates either.
Example 7. Chen, Knox, Arora et al. (2016) ^[7]^ considered clustered binary responses arising from multi-center studies. They modeled and simulated clustered binary responses from the following GLMM
(18)
where K denotes the number of study sites, n number of subjects within each site (cluster size), Xkj is a binary variable indicating treatment condition assigned to the i th subject within the k th site, and is the random effect accounting for correlations between responses yki from subjects with the k th site. They considered testing the
(19)
where pj = E(yki | xki = j) is the mean response yki= (or proportion of yki = 1) for the treatment group. The parameters β0 and β1 in (18) are related to p0 and p1 in (19) as
Thus given (19) and σ^2^λ, we can solve the above for β0 and β1 (numerically).
Chen, Knox, Arora et al. (2016) ^[7]^ considered σ^2^λ = 0, σ^2^λ = 0.1 and σ^2^λ = 1, with K = 20 clusters and n = 25 within each cluster. For each σ^2^λ, they obtained β0 and β1 and simulated data from the GLMM in (18) under Ha : p1 = 0.05. They fit the model (18) to the data generated using SAS NLMIXED, tested the null H0 : p0= 0.10 and rejected the null if the p-value is larger than α = 0.05. The process was repeated 1,000 times and power was estimated by the percent of times the null was rejected.
Shown in Table 4 are power estimates under “NLMIXED” for the three cases of σ^2^λ, along with true power values under (“True”) obtained using a different approach developed in Chen, Knox, Arora et al. (2016)^[7]^. As seen, power estimates from obtained from the SAS NLMIXED were quite different from the true power for each case of σ^2^λ, all underestimating the true power. Note that when σ^2^λ = 0, there is no data cluster and power only depends on the total sample size Kn. In this case, we can also obtain power estimates by using instead the SAS LOGISTIC NLMIXED for non-clustered data. This is indeed the case, since Chen, Knox, Arora et al. (2016) ^[7]^ reported that they obtained power estimates similar to 0.561 when running the Monte Carlo substitution to estimate power using the SAS LOGISTIC.
In this report, we discussed the two most popular regression models for longitudinal data. We focused on interpretation and computation of model parameters. For parameter interpretation, we discussed differences between the GLMM and WGEE when applied to model binary and count responses. Since parameters from the two longitudinal models are generally quite different, we should not expect similar estimates when applying the two models to real study data. Moreover, except for some special cases, it is generally not possible to find a relationship between estimates from the two models. Our analysis indicates that GLMM estimates can be quite difficult to interpret, while WGEE estimates afford straightforward interpretation.
Another major issue with the GLMM is to obtain reliable estimates using existing software. It seems that even software giants like SAS cannot provide correct estimates when applying GLMM to model binary responses. Until the computational problem is resolved, one may want to consider applying WGEE when modeling longitudinal binary responses.
We focused on the binary and count response when discussing interpretation of model parameters in this report. When modeling continuous responses, the two longitudinal models have the same interpretation for their model parameters and thus the interpretational issue does not arise. The computational problem seems only relevant to binary responses. For continuous responses, the log-likelihood function can be solved accurately. For count responses, major software such as R and SAS seem to provide quite reliable estimates.