Generalized estimating equations (GEEs)

McCullagh and Nelder (1989) introduced the Generalized Linear Model (GLM) for exponential family data with the form

fy (y, θ, φ) = exp｛f(yθ - b(θ))/a(φ) + c(y, φ)｝,

where a( ), b( ), and c( ) are given, θ is the canonical parameter, and φ is the dispersion parameter. The GLM is then given by

g(μ_i) = g(E[Y_i]) =x_i^'β,

where x_i is a p x 1 vector of covariates for the i^th subject, and β is a p x 1 vector of regression parameters. One of the attractive properties of the GLM is that it allows for linear as well as non-linear models under a single framework. It is possible to fit models where the underlying data are normal, inverse Gaussian, gamma, Poisson, binomial, geometric, and negative binomial by suitable choice of the link function g( ) (Hilbe, 1994).

Liang and Zeger (1986) and Zeger and Liang (1986) introduced generalized estimating equations (GEEs) to account for the correlation between observations in generalized linear regression models. One aspect of the approach builds upon previous methods of variance estimation developed to protect against inappropriate assumptions about the variance (Huber 1967; White 1980, 1982). The GEEs method is an extension of the generalized linear model, and is applicable in the analysis of correlated discrete outcome data (Liang and Zeger, 1986).

The appeal of GEEs is the interpretation of the data is the same as when a model that assumes independence is used, yet it is a valid method for correlated data (Zeger and Liang, 1992).

GEEs are often used when analyzing longitudinal or nested data. Several previous studies have examined the use of GEEs in the analysis of longitudinal injury-related and illness-related data. For example, Williamson et al. (1996) found when modeling longitudinal injury data under certain circumstances the use of GEEs resulted in different conclusions compared to logistic regression. While Williamson’s study addressed correlation of

in time. Diggle, Liang, and Zeger (1994) provided a detailed review of marginal models as well as other approaches, which including random effects models and transition models.

Let , i = 1,….., n; j = 1,……, t be the jth outcome for the ith subject, where we assume that observations on different subjects are independent, though we allow for association between outcomes observed on the same subject. In the GEEs setting, we are not assuming that is a member of the exponential family, but we are assuming that the mean and variance are characterized as in the GLM. We assume the marginal regression model

Yij

Yij

g(E[Y_ij]) = x_i^'β,

where x_ij is a p x 1 vector of study variables (covariates) for the i^th subject at the j^th outcome, β consists of the p regression parameters of interest and g( ) is the link function. Common choices for the link function might be g(a) = a for measured data (the identity link) g(a) = log(a) for count data (log link), or g(a) = log(a/(1-a)) for binary data (logit link). GEEs have been a popular approach to regression model fitting for this type of data. In the case of a binary data with the logit link, it will be that

log(E[Y_ij]/(1 - E[Y_ij])) = x_i^'β,

which implies that

E[Y_ij] = μ_ij= exp(x_ij^'β)/(1 + exp(x_ij^'β) ,

and if the outcomes are binary, it will be that var(Y_ij) =V_ij= exp(x_ij^' )/(1 + exp(x_ij^' )) ²,

In addition to this marginal mean model, it is needed to model the covariance structure of the correlated observations on a given subject. Assuming no missing data, the t x t covariance matrix of Y_i is modeled as

V = i φA_i^1/2R(α)A_i^1/2,

where A_i is a diagonal matrix of variance functions v(μ_ij), and R(α) is the working correlation matrix of Y _i indexed by a vector of parametersα.

2.3.2.2 Specification of working correlation matrix

There are a variety of common structures that may be appropriate to use to model the working correlation matrix. In general if the number of observations per cluster is small in a balanced and complete design, then an unstructured matrix is recommended. For datasets with mistimed measurements, it may be reasonable to consider a model where the correlation is a function of the time between observations. For datasets with clustered observations, there may be no logical ordering for observations within a cluster and an exchangeable structure may be most appropriate. Comparisons of estimates and standard errors from several different correlation structures may indicate sensitivity to misspecification of the variance structure.

For both the independence working structure and the fixed working structure, no estimation of α is performed. It is noted that use of the exchangeable working correlation matrix with measured data and identity link function is equivalent to a random effects model with a random intercept per cluster.

2.3.2.3 Empirical and model based variance estimators

Zeger and Liang (1986) referred to Vi as a working matrix because it is not required to be correctly specified for the parameter estimates and the estimated variance of the parameter estimates to be consistent (as long as the mean model itself is correct and there is no missing data). However, Liang and Zeger (1986) showed that there could be important gains in efficiency realized by correctly specifying the working correlation matrix. Sets of estimating equations are solved, through an iterative process, to find the value of the estimatorβ. An empirical variance estimator can be used to estimate var(β). This variance estimator is also referred to as a robust estimator. Another variance estimate available from GEE models is the

model are correctly specified. Since in general the analyst will not know the correct covariance structure, the empirical variance estimate will be preferred when the number of clusters is large. When the number of clusters is small, say < 20; the model based variance estimator may have better properties (Prentice 1988) even if the working variance is wrong.

This is because the robust variance estimator is asymptotically unbiased, but could be highly biased when the number of clusters is small.

2.3.2.4 Missing data issues

Longitudinal or clustered studies often have missing data, either by design or happenstance. If a litter in a teratology study is the level of clustering, litter size may vary between litters. If the missingness can be thought of as being missing completely at random in the sense of Little and Rubin (1987), then the consistency results established by Liang and Zeger (1986) hold. However, the notation and calculations for arbitrary missing data patterns are more complicated than in the balanced and complete case. Robins, Rotnitzky, and Zhao (1995) proposed methods to allow for data that is missing at random. Their inverse probability censoring weight approach requires that the missingness law be modeled, and that weights corresponding to the inverse probability of missingness be included in the GEEs. This will yield consistent parameter estimates, but the variance will tend to be incorrect.

Unfortunately, the method of Robins et al. (1995) only works well when there is dropout.

That is, once a subject misses a time, that subject is not seen again. Often subjects miss a single observation, and then are seen at the next time. In summary, when fitting GEEs, the analyst must consider not only the model for the mean, but also the model for the variance and the underlying missingness process.