The R Package geepack for Generalized Estimating Equations

1 JSSJ ournal of Statistical SoftwareJanuary 2006, Volume 15, Issue Generalized EstimatingEquationsUlrich HalekohDanish Institute of Agricultural SciencesS ren H jsgaardDanish Institute of Agricultural SciencesJun YanUniversity of IowaAbstractThis paper describes the core features of theRpackagegeepack, which implements thegeneralized Estimating Equations (GEE) approach for fitting marginal Generalized linearmodels to clustered data. Clustered data arise in many applications such as longitudinaldata and repeated measures. The GEE approach focuses on models for the mean of thecorrelated observations within clusters without fully specifying the joint distribution ofthe observations. It has been widely used in statistical practice. This paper illustrates theapplication of the GEE approach withgeepackthrough an example of clustered : Generalized Estimating equation, random effect, mixed model, IntroductionGeneralized Estimating Equations (GEE) (Liang and Zeger1986) are a general method foranalyzing data collected in clusters where 1) observations within a cluster may be correlated,2) observations in separate clusters are independent, 3) a monotone transformation of theexpectation is linearly related to the explanatory variables and 4) the variance is a functionof the expectation.

2 It is essential to note that the expectation and the variance referred to inpoints 3) and 4) are conditional given cluster-level or individual-level are several approaches to extend Generalized linear models to clustered data. Mixedeffect models and transition models (Diggle, Liang, and Zeger1994, Chapter 7, 9-10) fullyspecify the joint distribution within clusters via latent variables or conditional the presence of random effects, likelihood estimation necessitates the integration overthe random effects distributions, which may be numerically and Nelder(1996,2 TheRPackagegeepackfor Generalized Estimating Equations2001) introduced hierarchical Generalized linear models and showed that the integration maybe avoided by working on the h-likelihood. Compared to these approaches, the method ofGEE fits marginal mean models with the advantage that only correct specification of marginalmeans is needed for the parameter estimator to be consistent and asymptotically normal.

3 Thisapproach has become an important tool in analyzing longitudinal data or repeated measuresarising in a wide variety of applications. For a discussion on the relation between marginaland mixed effects models, seeHeagerty and Zeger(2000) andNelder and Lee(2004).Several implementations of GEE have become available (Horton and Lipsitz1999). Thebasic approach ofLiang and Zeger(1986) is available inSAS(SASI nstitute ,procgenmod),Stata(StataCorp LP2005),XLISP-STAT(Lumley1996) and inS-PLUSby thepackagesoswald(Smith1998) andgeeoryags(Carey2002,2004). The last two packageshave been ported toR(RDevelopment Core Team2005). TheRpackagegeepackimplementsthe basic approach and some extensions (Yan2002;Yan and Fine2004). Three features ofgeepackdistinguish it from other implementations: 1) There is an interface functiongeeglmwhich is designed to be as similar toglmas possible; 2) A jackknife variance estimator isavailable as an alternative to the sandwich estimator; and 3) Covariates can be incorporatedinto the scale and correlation parameters in a similar fashion to the mean modeling.

4 In thispaper, we illustrate the aspects ofgeepackwith the focus on the first two paper is organized as follows. In Section2we introduce an example dataset on repeatedmeasures of binary data. In Section3we outline the GEE approach and in Section4wedescribe the features of thegeeglmfunction that implements the approach ingeepack. Weclose the article with an analysis of the data in Section5and a conclusion in An example data setTo illustrate the type of problems well suited for the GEE approach we consider a data set onrespiratory illness. The data is provided ingeepackand detailed information about the datacan be found inKoch, Carr, Amara, Stokes, and Uryniak(1990). Briefly, the data comesfrom a clinical study in which the effect of a treatment of patients with respiratory illness wasVisitResponse pattern10 1 0 0 0 1 1 1 0 0 1 1 1 0 120 0 1 0 0 1 0 0 1 0 1 1 0 1 130 0 0 1 0 0 1 0 1 1 1 0 1 1 140 0 0 0 1 0 0 1 0 1 0 1 1 1 1 Baseline TreatSum0A7 2 2 2 1 0 1 0 1 0 1 2 0 4 7 30P 18 1 0 2 1 2 0 0 1 0 0 1 2 0 3 311A0 0 0 0 0 0 1 1 0 0 4 0 1 0 17 24P1 4 1 0 0 0 0 1 1 3 1 1 2 1 10 26 Sum26 7 3 4 2 2 2 2 3 3 6 4 5 5 37 111 Table 1: Number of patients for the different response patterns classified by baseline re-sponse and treatment.

5 The patterns are ordered according to increasing numbers of of Statistical Software3 Number of positive reponses0 1 2 34 Sex F7 3 3 37M19 13 9 1730 Center 118 9 6 111228 7 6 925 Table 2: Number of patients for the number of positive responses across the four visits forsex and A total of 111 patients from two clinical centers were randomized to receive eitherplacebo or an active treatment. At four examination visits, the respiratory state of a patientwas classified as good (=1) or poor (=0). Variables characterizing a patient were:center(1,2), treatmenttreat(A=active, P=placebo),sex(M=male, F=female),age(in years) atbaseline andbaselinerespiratory state which is also binary. The values of the covariateswere constant for the repeated elementary observations on each the number of patients for the response patterns across the four visits classifiedby baseline status and treatment.

6 Patients with baseline respiratory status equal to 0 seemeither to have a low or large number of positive responses whereas patients with a baselineof 1 tend to respond positively. Table2describes the distribution of the number of positiveresponses per patient for sex and the plot of age against the proportion of positive responses for each indicates a quadratic relationship between the proportions and the age. We fit a logisticmodel to the data (which would be appropriate if there were no time effect and no spread inthe response probabilities for patients with the same covariate values):R> <- glm(outcome ~ baseline + center + sex + treat + age + I(age^2),+ data = respiratory, family = binomial)The correlation matrix of the Pearson residuals within a patient based on theglm-fit is shownin Table3and indicates an appreciable correlation within patient visit2 visit3 visit4visit1 3.

7 Correlation matrix for the measurements at different visits based on the Pearsonresiduals from the logistic Theory of GEEFor the regression analysis of correlated observations,Liang and Zeger(1986) introducedthe GEE approach. This approach Generalized the estimation method of quasi-likelihood of4 TheRPackagegeepackfor Generalized Estimating Equationsl llllllllllll lllllllll lll lll lllllll lllllllll llllllllllllllllllllllllllll lllllllllll lllllllllllllllllllll of positive responses per patient against 1: Relation ofageto the proportion of positive responses. The smooth line is a (1974) to correlated data. An alternative generalization was proposed byLeeand Nelder(1996,2001). An extensive review of the development of the GEE approach isgiven byZiegler, Kastner, and Blettner(1998).

8 Consider a sample ofi= 1, .. , Kindependent multivariate observationsYi= (Yi1, .. , Yit, .. , Yini). Hereimay represent a cluster withniobservations. The ex-pectationsE(Yit) = itare related to thepdimensional regressor vectorxitby the mean-linkfunctiongg( it) =x>it .(1)LetVAR(Yit) = ait,(2)where is a common scale parameter andait=a( it) is a known variance function. LetRi( ) be a working correlation matrix completely described by the parameter vector oflengthm. LetVi= A1/2iRi( )A1/2ibe the corresponding working covariance matrix ofYi, whereAiis the diagonal matrix withentriesait. For given estimates ( , ) of ( , ) the estimate is the solution of the equationK i=1 >i V 1i(Yi i) = and Zeger(1986) suggest to use consistent moment estimates for and . This yieldsan iterative scheme which switches between Estimating for fixed values of and andestimating ( , ) for fixed values of.

9 This scheme yields a consistent estimate for . More-over,K1/2( ) is asymptotically multivariate normally distributed with zero mean andJournal of Statistical Software5namefunctionv( )gaussianidentitybinomial (1 ), (0,1)poisson , >0gamma 2, >0 Table 4: Variance functions matrix = limK K 10 1 10where 0=K I=1 >i V 1i i >, 1=K i=1 >i V 1iCOV(Yi)V 1i i >.(3)Replacing , and by consistent estimates and the covariance matrixCOV(Yi) by (Yi i)(Yi i)>in (3) yields a so called sandwich estimate of . The estimate is a consistentestimate of even if the working correlation matricesRi( ) are ThegeeglmfunctionThegeeglmfunction largely follows the syntax of theglmfunction and many of the methodsavailable forglmobjects are also available forgeeglmobjects.

10 We discuss in the followingthe most important arguments of thegeeglmfunction and theanovamethod for comparingmodels by Wald tests. These will be exemplified in Variance and link functions (family)The variance function is specified by thefamilyargument ingeeglmand is identified by thename of the corresponding distribution in a Generalized linear model. The available variancefunctions are given in Table4. The available link functions for the mean are the same asthose inglmwith the exception of thecauchitlink for Working correlation (corstr, zcor)Four pre-defined working correlation structures are available and are specified via the argu-mentcorstr(Table5). It is also possible to provide a correlation matrix the entries of whichremain fixed under the computation.