Generalized Estimating Equations (gee) for glm–type data

1 Generalized Estimating Equations (gee)for glm type dataS ren H Research UnitDanish Institute of Agricultural SciencesJanuary 23, 2006 Printed: January 23, 2006 File: Preliminaries32 Working example respiratory illness43 Correlated Pearson residuals94 Marginal vs. conditional models125 Marginal models for glm type data146 Estimating Equations for gee type Specifications needed for Deriving and solving Newton Estimation of the covariance of .. Model based Emperical estimate sandwich The working correlation Exploring different working correlations338 Comparison of the parameter When do GEEs work?.. What to do geeglm vs. 23, 2006page 231 PreliminariesThese notes deal with fitting models for responses of type oftendealt with with Generalized linear models (glm) but with thecomplicating aspect that there may be repeated measurementson the same approach here is Generalized Estimating Equations (gee).There are two packages for this purpose in R: geepack and focus on the former and note in passing that the latter doesnot seem to undergo any further geepack package is described in the paper by Halekoh,H jsgaard and Yun in Journal of Statistical Software, , January 23, 2006page 342 Working example respiratory illnessExample 1 The data are from a clinical trial of patients withrespiratory illness, where 111 patients from two different clinicswere randomized to receive either placebo or an activetreatment.

2 Patients were examined at baseline and at four visitsduring treatment. At each examination, respiratory status(categorized as 1 = good, 0 = poor) was determined. The recorded variables are:Center (1,2), ID, Treatment (A=Active, P=Placebo),Gender (M=Male,F=Female), Age (in years at baseline),Baseline Response. The response variables are:Visit 1 Response, Visit 2 Response, Visit 3 Response, Visit4 for 8 patients are shown in Table 23, 2006page 45center id treat sex age baseline visit1 visit2 visit3 visit411 1 PM460000021 2 PM280000031 3 AM231111141 4 PM441111051 5 PF131111161 6 AM340000071 7 PM430101181 8 AM2800000 Table 1: Respiratory data for eight individuals. Measurements onthe same individual tend to be 23, 2006page 56 Interest is in comparing the treatments, but also to includecenter, age, gender and baseline response in the Table 1 it is clear, that there is a dependency among theresponse measurements on the same person measurements onthe same person tend to be dependency must be accounted for in the modelling.

3 January 23, 2006page 67 Example 2A first approach is to ignore the dependency. Thisapproach isnotappropriate but the response measured on theith person at visitv, wherev= 1,..,4 Since the response outcomes are binary,yiv {0,1}, it is tempting to consider the binomial distribution asbasis for the modelling. That is, to assume thatyiv bin(1, iv)and that allyivare specification of ivwe consider in the following the linearpredictorlogit( iv) = + center(i)+ treat(i)+ agei+ baselineiNote that the expression forlogit( iv)does not include the vistv. We will write this briefly aslogit( ) =center+treat+age+baselineOther linear predictors can clearly be 23, 2006page 78 Table 2 contains the parameter estimates under the Std. Error z value Pr(>|z|)(Intercept) 2: Parameter estimates when assuming independenceA more elaborate model is to allow for is, to make a quasi likelihood model where the variancefunction is iv(1 iv)The dispersion parameter is estimated as = , thereis no indication of over/under dispersion.

4 January 23, 2006page 893 Correlated Pearson residualsBased on the fitted (wrong) independence model we cancalculate the Pearson residualseiv=yiv ivp iv(1 iv), i= 1,..,N,v= 1,..,4which under the model approximately have mean 0 and these we can estimate the covariance matrix in Table 3which shows covariances between measurements on the 23, 2006page 91012341 3: Covariance matrix based on Pearson 23, 2006page 1011 Since the elements on the diagonal in Table 3 are about 1, thematrix can also be regarded as acorrelation the observationswereindependent then the true ( ) correlations should be estimated correlations in Table 3 suggest that there is apositive task in the following is to account for the correlationbetween measurements on the same 23, 2006page 11124 Marginal vs. conditional modelsLinear mixed models of the typey=X +Zu+ecan be described as conditional models. Suppose thatVar(e) = 2I. Then the conditional distribution ofygivenuisy|u N(X +Zu, 2I)so in the conditional distribution, the components , we impose a structure onuin terms ofVar(u) = marginal distribution ofyisy N(X ,ZGZ+ 2I)so marginally the components ofyare dependent with thestructure given inV=ZGZ + 23, 2006page 1213So this way, one can see the linear mixed model formula as a wayof building up a model in which the responses are alternative approach is to construct a marginal modeldirectly, N(X ,V)by specifying directly a structure 23, 2006page 13145 Marginal models for glm type dataA way of dealing with correlated glm type observations is tocreate a marginal model directly.

5 That is, to create a modelfor a 4 dimensional response vector which consist of follow the approach by Liang and Zeger (1986).The term marginal model is quoted, because formally we donot specify a proper statistical model (in terms of makingdistributional assumptions).All we do is to to the meanE(y)depends on the covariates through a linkfunctiong(E(y)) =X the varianceVar(y)varies as a function of the mean(the variance function), (y) =v(E(y)).January 23, 2006page 1415 With these specifications, one can derive a system of estimatingequations by which an estimate andVar( )can be N( ,..).So we make fewer assumptions than if we specify a fullstatistical model. This extra flexibility comes at a price: The estimate may not be the best possible. Hypothesis testing is based on Wald tests (since, as there isno distribution and hence no likelihood). Model checking is can regard GEEs as a quick and dirty 23, 2006page 15166 Estimating Equations for gee typedataFor correlated glm type data, Estimating Equations have in thelitterature become known asgeneralised Estimating Equations (GEEs).

6 GEEs can, in connection with correlated glm type data, beregarded as an extension of the esimation methods (scoreequations) used GLMs/QLs. This justifies the term Generalized . On the other hand, the Estimating Equations used inconnection with correlated glm type data are are ratherspecialized type of Estimating Equations . As such, the term Generalized is a little this reason the function for dealing with these types of datain the geepack package is called geeglm().January 23, 2006page 1617 With GEEs for GLM type data the emphasis is on modeling the expectation of thedependent variable in relation to the covariates (just likewith GLMs), whereas the correlation structure is considered to be anuisance(notof interest in itself), which is accounted for by the 23, 2006page Specifications needed for GEEsThe setting is as follows: On each ofi= 1,..,Nsubjects, thereare madenimeasurementsyi= (yi1,..,yini). Measurements on different subjects are assumed to beindependent Measurements on the same subject are allowed to model formulation is similar to that of a GLM:Systematic part:Relate theexpectationE(yit) = itto thelinear predictor via the link functionh( it) = it=x it Random part:Specify how the varianceVar(yit)is related tothe meanE(yit)by specifying avariance functionV( it)such thatVar(yit) = V( it).

7 January 23, 2006page 1819 The correlation part:In addition to these GLM steps weneed to impose a correlation structure for observations onthe same unit. This is done by specifying aworking correlation 23, 2006page Deriving and solving GEEsConsider observationsy1,..,ynwith common mean . The leastsquares criterion for Estimating is to minimize ( ,y) =Xi(yi )2 This is achieved by setting the derivative to zero: ( ;y) = ( ;y) = 2Xi(yi ) = 0We say that ( ) = ( ;y)is an Estimating function and ( ;y) = 0is an Estimating 23, 2006page 2021 Consider datay= (y1,..,yn)and a modelp(y; )where Rp. The idea behind Estimating functions is to find afunction ( )which immitates the score functionU( ) =dd logp(y; ).Let ( )= ( ,y)be such a function denoted an Estimating function. Solve(usually by iteration) the Estimating Equations ( ) = 0giving = (y)IfE ( ( )) = 0for all (which holds for the score function), then is said to be unbiased. (Note that unbiasedness is a propertyof the estimation function rather than of the estimate.)

8 January 23, 2006page 2122In practice we frequently we consider weighted sums ofestimating functions of the formXiai(yi i( ))which consequently is 23, 2006page 2223 Unbiasedness of the estimation function implies that isasymptotically consistent. Another property is that isasymptotically normal. The sensitivity of the Estimating function is thep pmatrixS ( ) =E ( )which indicates how steep is on average . So a largevalue ofSis good. The variability of an Estimating function isV ( ) =V ar ( ( )) =E ( ( ) ( )>)A small value ofVis good because that indicates thatdifferent samples give almost the same . The Godambe information matrix is defined asJ ( ) =S> ( )V ( ) 1S ( )January 23, 2006page 2324 It then holds that approxN( ,J ( ) 1)Note: If is the score function (arising from a likelihood) thenS ( ) = I( )andV ( ) =I( )and henceJ ( ) =I( ).If ( )has the form =X>(y )thenS ( ) =X>[ 1: : p].Moreover,V ( ) =E( ( ))2=E(X>(y )(y )>X)which is not so easy to calculate.

9 In principle, however, thisquantity can be estimated using V ( ) =X>(y )(y )>XIn practice, this may cause problem since V in this case maynot be 23, 2006page 2425 The GEE by Liang and Zeeger (1986) for Estimating apvector is given by ( ) =Xi i 1i(yi i( ))(1)Sinceg( ij) =x ij , the deriavative i has as itskjth entry[ i ]kj=xikjg ( ik)The variance is i= A1/2iR( )A1/2iwhereAiis diagonal with thev( ij)s on the diagonal andR( )isthe correlation correlation matrix is generally unknown, so therefore onespecifies a working correlation matrix , with anautoregressive in a repeated measurements absence of a good guess ofR( )the identity matrix is often agood 23, 2006page 2526 TypicallyR( )is estimated from data iteratively by using thecurrent estimate of to calculate a function of the Pearsonresidualseij=yij ij( )pv( ij( ))(2)The dispersion parameter is often estimated as =1N pNXi=1niXj=1e2ij(3)January 23, 2006page Newton iterationThe fitting algorithm then an initial estimate of from a glm ( byassuming independence) an estimateR( )of the working correlation onthe basis of the current Pearson residuals and the currentestimate of an estimate of the variance as i= A1/2i R( )A1 an updated estimate of based on theNewton step := + [Xi i 1i i ] 1[Xi i 1i(yi i( ))]January 23, 2006page 2728 Iterate through 2 4 until convergence.

10 Note that needs not tobe estimated until the last GEE estimate of is often very similar to the estimateobtained if observations were treated as being independent. Inother words, the estimate for GEEs is often very similar to theestimate obtained by fitting a QL model to the 23, 2006page Estimation of the covariance of There are two classical ways of Estimating the covarianceCov( ). Model based estimateCov( )m=I 10, I0=Xi i 1i i This is the GEE version of inverse Fisher information usedwhen esitmatingCov( )in a glm. HereCov( )mconsistentlyestimatesCov( )if i) the mean model and ii) the workingcorrelation are 23, 2006page Emperical estimate sandwich estimateCov( )e=I 10I1I 10whereI1=Xi i 1iCov(yi) 1i i HereCov( )eis a consistent estimate ofCov( )even if theworking correlation is misspecified, ifCov(yi)6= practice,Cov(yi)is replaced by(yi i( ))(yi i( )) .January 23, 2006page The working correlation matrixThe notion of working correlation matrices can be introducedthrough Example 1 on the respiratory that the measurements on theith person areyi= (yi1.