Generalized Estimating Equations
Let Y_{ij}, j = 1, ... ,n_{i}, i = 1, ... ,K
represent the jth measurement on the ith subject.
There are n_{i} measurements on subject i
and total measurements.
Correlated data are modeled using the same link
function and linear predictor setup (systematic
component) as the independence case.
The random component is described by the same variance
functions as in the independence case, but the covariance
structure of the correlated measurements must also be modeled.
Let the vector of measurements on the ith subject be
Y = [Y_{i1}, ... , Y_{i ni}]' with
corresponding vector of means and let V be
the covariance matrix of Y.
Let the vector of independent, or explanatory, variables
for the jth measurement on the ith subject be

X_{ij} = [x_{ij1}, ... , x_{ijp}]'
The Generalized Estimating Equation of Liang and Zeger (1986)
for estimating the p ×1 vector of regression parameters
is an extension of the independence
estimating equation to correlated data and is given by
Since
where g is the link function, the p ×n_{i} matrix
of partial derivatives of the mean with respect to the
regression parameters for the ith subject is given by
Working Correlation Matrix
Let be an n_{i} ×n_{i}
"working" correlation matrix that is fully
specified by the vector of parameters .The covariance matrix of Y is modeled as
where A_{i} is an n_{i} ×n_{i} diagonal matrix
with as the jth diagonal element.
If is the true correlation matrix of
Y, then V is the true covariance matrix of Y.
The working correlation matrix is usually unknown and must be estimated.
It is estimated in the iterative fitting process using
the current value of the parameter vector to
compute appropriate functions of the Pearson residual
If you specify the working correlation as
R_{0} = I, which is the identity matrix, the
GEE reduces to the independence estimating equation.
Following are the structures of the working
correlation supported by the GENMOD procedure and the
estimators used to estimate the working correlations.
Working Correlation Structure


Estimator

Fixed   
 Corr(Y_{ij},Y_{ik}) = r_{jk}
where r_{jk} is the jkth element of
a constant, userspecified correlation matrix
R_{0}.  The working correlation is not estimated
in this case. 
  
Independent   
  The working correlation is not estimated
in this case. 
  
mdependent   
  
  
  
Exchangeable   
  
  
  
Unstructured   
  
  
Autoregressive AR(1)   
 for t = 0,1,2, ... ,n_{i}j  
  
Dispersion Parameter
The dispersion parameter is estimated by
where is the total number of measurements
and p is the number of regression parameters.
The square root of is reported
by PROC GENMOD as the scale parameter in the
"Analysis of GEE Parameter Estimates
ModelBased Standard Error Estimates" output table.
Fitting Algorithm
The following is an algorithm for
fitting the specified model using GEEs.
Note that this is not in general a likelihoodbased method
of estimation, so that inferences based on likelihoods are
not possible for GEE methods.
 Compute an initial estimate of with
an ordinary generalized linear model assuming independence.
 Compute the working correlations R
based on the standardized residuals, the current , and
the assumed structure of R.
 Compute an estimate of the covariance:
 Update :
 Iterate steps 24 until convergence
Missing Data
Refer to Diggle, Liang, and Zeger (1994, Chapter 11) for a
discussion of missing values in longitudinal data.
Suppose that you intend to take measurements
Y_{i1}, ... , Y_{in} for the ith unit.
Missing values for which Y_{ij} are missing whenever Y_{ik}
is missing for all are called dropouts.
Otherwise, missing values that occur intermixed with
nonmissing values are intermittent missing values.
The GENMOD procedure can estimate the working correlation
from data containing
both types of missing values using the all
available pairs method, in which all nonmissing pairs
of data are used in the moment estimators of the
working correlation parameters defined previously.
For example, for the unstructured working correlation model,
where the sum is over the units that have nonmissing
measurements at times j and k, and K' is the
number of units with nonmissing measurements at j and k.
Estimates of the parameters for other working correlation
types are computed in a similar manner, using available
nonmissing pairs in the appropriate moment estimators.
The contribution of the ith unit to the parameter
update equation is computed by omitting the
elements of , the columns of
, and the rows
and columns of V corresponding to missing measurements.
Parameter Estimate Covariances
The modelbased estimator of
is given by
where
This is the GEE equivalent of the inverse of the Fisher
information matrix that is often used in generalized
linear models as an estimator of the covariance estimate
of the maximum likelihood estimator of .It is a consistent estimator of the covariance
matrix of if the mean model and the
working correlation matrix are correctly specified.
The estimator
is called the empirical, or robust, estimator
of the covariance matrix of where
It has the property of being a consistent estimator
of the covariance matrix of , even if
the working correlation matrix is misspecified, that
is, if .In computing M, and are replaced by
estimates, and Cov(Y) is replaced by the estimate
Multinomial GEEs
Lipsitz, Kim, and Zhao (1994) and Miller, Davis, and Landis (1993)
describe how to extend GEEs
to multinomial data. Currently,
only the independent working correlation is available
for multinomial models in PROC GENMOD.
If the responses are binary (that is, they take only two values),
then there is an alternative method to account for the association
among the measurements. The Alternating Logistic Regressions (ALR)
algorithm of Carey, Zeger, and Diggle (1993) models the association
between pairs of responses with log odds ratios, instead of with
correlations, as ordinary GEEs do.
For binary data, the correlation between the jth and kth
response is, by definition,
The joint probability in the numerator satisfies the following
bounds, by elementary properties of probability, since
:
The correlation, therefore, is constrained to be within limits that
depend in a complicated way on the means of the data.
The odds ratio, defined as

OR(Y_{ij},Y_{ik}) = [(Pr(Y_{ij} = 1,Y_{ik} = 1)Pr(Y_{ij} = 0,Y_{ik} = 0))/(Pr(Y_{ij} = 1,Y_{ik} = 0)Pr(Y_{ij} = 0,Y_{ik} = 1))]
is not constrained by the means and is preferred,
in some cases, to correlations for binary data.
The ALR algorithm seeks to model the logarithm of the odds ratio,
, as
where is a q ×1 vector of regression parameters
and z_{ijk} is a fixed, specified vector of coefficients.
The parameter can take any value in with
corresponding to no association.
The log odds ratio, when modeled in this way with a regression model,
can take different values in subgroups defined by z_{ijk}.
For example, z_{ijk} can define subgroups within clusters, or
it can define `block effects' between clusters.
You specify a GEE model for binary data using log odds ratios
by specifying a model for the mean, as in ordinary GEEs, and
a model for the log odds ratios. You can use any of the
link functions appropriate for binary data in the model for the
mean, such as logistic, probit, or complementary loglog.
The ALR algorithm alternates between a GEE step to update the
model
for the mean and a logistic regression step to update the log odds ratio
model.
Upon convergence, the ALR algorithm
provides estimates of the regression parameters
for the mean, , the regression parameters for the log odds ratios,
, their standard errors, and their covariances.
Specifying Log Odds Ratio Models
Specifying a regression model for the log odds ratio requires
you to specify rows of the zmatrix
z_{ijk} for each cluster i and each unique withincluster
pair (j, k).
The GENMOD procedure provides several methods of specifying z_{ijk}.
These are controlled by the LOGOR=keyword and associated
options in the REPEATED statement.
The supported keywords and the resulting log odds ratio models
are described as follows.
 EXCH
 exchangeable log odds ratios. In this model, the log odds ratio
is a constant for all clusters i and pairs (j,k). The parameter
is the common log odds ratio.

z_{ijk} = 1 for all i, j, k
 FULLCLUST
 fully parameterized clusters.
Each cluster is parameterized in the same way, and
there is a parameter
for each unique pair within clusters. If a complete cluster is of size
n, then there are [(n(n1))/2] parameters in the vector
. For example, if a full cluster is of size 4,
then there are [4 ×3/2] = 6 parameters,
and the zmatrix is of the form
The elements of correspond to log odds ratios for cluster
pairs in the following order:
Pair

Parameter

(1,2)  Alpha1 
(1,3)  Alpha2 
(1,4)  Alpha3 
(2.3)  Alpha4 
(2,4)  Alpha5 
(3,4)  Alpha6 
 LOGORVAR(variable)
 log odds ratios by cluster.
The argument variable is a variable name
that defines the `block effects' between clusters. The log odds ratios
are constant within clusters, but they take a different value for
each different value of the variable. For example, if
Center is a variable in the input data set taking a different
value for k treatment centers, then specifying
LOGOR=LOGORVAR(Center) requests a model with different
log odds ratios for each of the k centers, constant within center.
 NESTK
 knested log odds ratios. You must also specify the
SUBCLUST=variable option
to define
subclusters within clusters. Within each cluster, PROC GENMOD computes
a log odds ratio parameter for pairs having the same value
of variable for both members of the
pair and one log odds ratio parameter for each unique combination
of different values of variable.
 NEST1
 1nested log odds ratios. You must also specify the
SUBCLUST=variable option to define
subclusters within clusters. There are
two log odds ratio parameters for this model. Pairs having the same
value of variable correspond to one parameter;
pairs having different
values of variable correspond to the other parameter.
For example, if clusters are hospitals and subclusters are
wards within hospitals, then patients within the same ward
have one log odds ratio parameter, and patients from different wards
have the other parameter.
 ZFULL
 specifies the full zmatrix. You must also specify
a SAS data set containing the zmatrix with
the ZDATA=datasetname option. Each observation in the data set
corresponds to one row of the zmatrix. You must specify the ZDATA
data set as if all clusters are complete, that is, as if all
clusters are the same size and there are no missing observations.
The ZDATA data set has
K [n_{max}(n_{max}1)/2] observations, where K is the number of
clusters and n_{max} is the maximum cluster size.
If the members of cluster i are ordered as 1,2, ... ,n, then
the rows of the zmatrix must be specified for pairs in the order
(1,2),(1,3), ... ,(1,n),(2,3), ... ,(2,n), ... ,(n1,n).
The variables specified in the REPEATED statement
for the SUBJECT effect must also be present in the ZDATA=
data set to identify clusters.
You must specify variables in the data set
that define the columns of the zmatrix by the
ZROW=variablelist option. If there are q columns,
(q variables in variablelist), then there are q
log odds ratio parameters. You can optionally specify
variables indicating the cluster pairs corresponding to each
row of the zmatrix with the YPAIR=(variable1, variable2)
option.
If you specify this option, the data from the ZDATA data set is
sorted within each cluster by variable1 and
variable2.
See Example 29.6 for an example of specifying a full
zmatrix.
 ZREP
 replicated zmatrix. You specify zmatrix data exactly
as you do for the ZFULL case, except that you specify only one complete cluster.
The zmatrix for the one cluster is replicated for each cluster.
The number of observations in the ZDATA data set is
[(n_{max}(n_{max}1))/2], where n_{max} is the size of a complete
cluster (a cluster with no missing observations).
See Example 29.6 for an example of specifying a replicated
zmatrix.
 ZREP(matrix)
 direct input of the replicated zmatrix.
You specify the zmatrix
for one cluster with the syntax
LOGOR=ZREP ( (y_{1} y_{2}) z_{1} z_{2} ... z_{q}, ... ), where y_{1} and y_{2}
are numbers representing a pair of observations and
the values z_{1}, z_{2}, ... , z_{q}
make up the corresponding row of the zmatrix.
The number of rows specified is
[(n_{max}(n_{max}1))/2], where n_{max} is the size of a complete
cluster (a cluster with no missing observations).
For example,
LOGOR = ZREP((1 2) 1 0,
(1 3) 1 0,
(1 4) 1 0,
(2 3) 1 1,
(2 4) 1 1,
(3 4) 1 1)
specifies the [4 ×3/2] = 6
rows of the zmatrix for a cluster of size 4 with q=2
log odds ratio parameters. The log odds ratio for pairs
(1 2), (1 3), (1 4) is , and the log odds ratio for
pairs (2 3), (2 4), (3 4) is .
Boos (1992) and Rotnitzky and Jewell (1990) describe score
tests applicable to testing in GEEs,
where L is a userspecified r ×p contrast matrix or a
contrast for a Type 3 test of hypothesis.
Let be the regression parameters resulting
from solving the GEE under the restricted model
, and let be the
generalized estimating equation values at .
The generalized score statistic is
where is the modelbased covariance estimate and
is the empirical covariance estimate.
The pvalues for T are computed based on the chisquare distribution
with r degrees of freedom.
Copyright © 1999 by SAS Institute Inc., Cary, NC, USA. All rights reserved.