Fits models to longitudinal data by generalized estimating equations (D.M. Smith & M.G.Kenward).
Options
PRINT = string token |
What to display (estimates, correlations, scalefactor, wald, monitoring); default esti, corr, scal |
|---|---|
DISTRIBUTION = string token |
Distribution of response (normal, Poisson, binomial, gamma, inversenormal, negativebinomial); default * |
LINK = string token |
Link function (identity, logarithm, logit, reciprocal, power, squareroot, probit, complementaryloglog, logratio); default * |
EXPONENT = scalar |
Exponent for power link; default -2 |
TERMS = formula |
Explanatory variates, factors etc |
CONSTANT = string token |
How to treat constant (estimate, omit); default esti |
FACTORIAL = scalar |
Limit for expansion of model terms; default 3 |
AGGREGATION = scalar |
Fixed parameter for negative binomial distribution (parameter k as in variance function var = mean + mean2/k); default 1 |
KLOGRATIO = scalar |
Parameter for logratio link, in form log(mean / (mean + k)); default as set in AGGREGATION option |
QUADESTIMATION = string token |
Whether to use quadratic estimation (used, notused); default used |
SCALEFACTOR = string token |
How to calculate the scale factor (fixed, constant, varytime); default varies with distribution, fixed for Poisson, binomial and negative binomial, constant for rest |
SFVALUE = scalar |
Value for scale factor when SCALEFACTOR=fixed; default 1.0 for Poisson and binomial, missing for rest |
CRTYPE = string token |
Form of correlation matrix (independence, unstructured, exchangeable, autoregressive, dependence, antedependence); default * |
ORDER = scalar |
Order in dependence and ante-dependence form of correlation matrix; default 1 |
TIMEDEPENDENT = string token |
Whether correlation in dependence model changes with time (no, yes); default no |
Parameters
Y = variates |
Response variate for each analysis |
|---|---|
NBINOMIAL = variates or scalars |
Denominator in binomial |
FITTEDVALUES = variates |
To store fitted values |
RESIDUALS = variates |
To store residuals |
SUBJECT = factors |
Identifier of subjects |
OUTCOME = factors |
Identifier of outcomes |
COUNT = variates |
Variate of counts of no. outcomes |
TIME = factors |
Times of repeated measures variate |
WEIGHT = variates |
Weight variate |
OFFSET = variates |
Offset variate |
SAVE = pointers |
Structure to save output variables |
Description
GEE implements the General Estimating Equation (GEE) methodology of Liang & Zeger (1986) with quadratic estimation for the covariance structure. In the terminology of Liang et al. (1992) the methodology implemented is a form of GEE1. Full details of the implementation are given in Kenward & Smith (1995a). GEE, as implemented here, is a comparatively simple non-likelihood method for fitting marginal models to repeated measurements that can be used when the response has a distribution in the exponential family. This includes the Gaussian distribution, for which the procedure implemented here reduces to a form of the EM algorithm, and then produces exact ML or REML estimates, or a close approximation to these depending on the particular correlation structure chosen. For other distributions the resulting estimates are not maximum likelihood but can be shown to have asymptotic properties familiar from quasi-likelihood, such as consistency and asymptotic normality.
The standard range of generalized linear models (as in procedure GLM) can be fitted involving a variety of covariance or correlation structures over the times of the repeated measurements. The standard links and distributions can be chosen by setting the options DISTRIBUTION, LINK, EXPONENT, AGGREGATION and KLOGRATIO, as in the MODEL directive. Non-standard ones require the definition of auxiliary procedures to carry out the necessary calculations (see the Method section). The terms in the fitted model are specified by the TERMS option, which may be set to a formula or left unset to fit a null model. The FACTORIAL option (default 3) sets a limit on the number of factors and variates in the terms that are fitted, as in the FIT directive. The CONSTANT option can be used to omit a constant term. Setting the QUADESTIMATION option to used requests the use of quadratic estimation for the data-based covariance or correlation matrix (see Kenward & Smith 1995a). The SCALEFACTOR option specifies the form of scalefactor to be used (fixed to a value specified by the SFVALUE option, constant over times of repeated measurements, or varying over times of repeated measurements). The CRTYPE option specifies the structure of the covariance or correlation matrix over the times of the repeated measurements. The ORDER option specifies the order of the covariance or correlation structures for the dependence and ante-dependence cases, with option TIMEDEPENDENT specifying whether the correlation in a dependence structure changes with the time of the repeated measurement.
The Y parameter must be set to specify the response variate. For a binomial distribution the NBINOMIAL parameter must also be set; this may be either a variate or a scalar (if it is a scalar, GEE maps it out automatically to a variate with the same number of values as Y). The SUBJECT parameter specifies a factor to identify the subjects. Alternatively, where the data consist of outcomes and numbers with those outcomes, the parameter OUTCOME must be set to the identifier of the outcome and the parameter COUNT to the number with the outcome. The parameter TIME must be set to the times of the repeated measurements. The parameters WEIGHT and OFFSET specify weight and offset variates that may be involved.
The output from the procedure is controlled by the PRINT option; by default estimates, their standard errors, covariances or correlations and scalefactors are given. Two sets of standard errors are provided for the estimates. One is the naive estimate which assumes the specified covariance or correlation structure holds. The other is the sandwich estimate which makes no such assumption. When PRINT=wald, Wald tests are produced using both sets of standard errors and correlations.
The fitted values and residuals can be obtained by setting the parameters FITTEDVALUES and RESIDUALS. The residuals are the Pearson residuals as defined in the Guide to the Genstat Command Language, Part 2, Section 3.1.1.
The SAVE parameter can save various details of the analysis, in a pointer with the following suffixes and labels:
1 or 'scalefactors' |
scalefactor(s), |
|---|---|
2 or 'correlation' or 'covariance' |
correlations or covariances, according to the type of model (and labelled appropriately), |
3 or 'estimates' |
the estimates of the linear predictor parameters, |
4 or 'naive covariances' |
naive variance-covariance matrix for the estimates, |
5 or 'sandwich covariances' |
sandwich variance-covariance matrix for the estimates, |
6 or 'naive Wald' |
Wald tests calculated using the naive variance-covariance matrix, and |
7 or 'sandwich Wald' |
Wald tests calculated using the sandwich variance-covariance matrix. |
The algorithms in the procedure have been set up assuming that the data contain a complete set of observations for each subject. Where there are missing values these must be included explicitly (using the missing value symbol *) to create a complete set of observations. Missing values are allowed in both the Y variate and the explanatory variates in TERMS.
In the case of the Gaussian distribution, a working covariance matrix, rather than correlation matrix, is used. This provides considerable simplification within the algorithm.
This is a complicated algorithm and some examples may take a while to run. If necessary, however, you can set option PRINT=monitoring to see what is happening.
Options: PRINT, DISTRIBUTION, LINK, EXPONENT, TERMS, CONSTANT, FACTORIAL, AGGREGATION, KLOGRATIO, QUADESTIMATION, SCALEFACTOR, SFVALUE, CRTYPE, ORDER, TIMEDEPENDENT.
Parameters: Y, NBINOMIAL, FITTEDVALUES, RESIDUALS, SUBJECT, OUTCOME, COUNT, TIME, WEIGHT, OFFSET, SAVE.
Method
For full details of the method implemented in this procedure see Kenward & Smith (1995a). A generalized linear model is formulated for the marginal distribution of the observations at each time point using an appropriate link function and error distribution. If the repeated measurements could be assumed to be independent, the well-known iterative weighted least squares fitting procedure could be used to obtain ML estimates of the marginal model parameters. However this ignores the dependence among the repeated measurements. Full likelihood is in general very awkward in this setting so, to avoid a formal introduction of dependence into the model, a working correlation matrix is introduced into the iterative procedure, changing the least squares from a weighted to a generalized form. The correlation matrix can be introduced in various ways. It can be held constant throughout the iterative procedure. An example of this is the use of the identity matrix, leading to the so-called independence estimating equations for which the process reduces back to that of fitting a univariate generalized linear model. Alternatively an estimated correlation matrix can be introduced into the algorithm which is updated at each cycle using quadratic estimation: essentially the correlation structure is estimated from the residuals using the equations that would be appropriate were the residuals normally distributed. On convergence consistent estimates of the marginal linear model parameters are obtained and, if the correlation structure chosen is appropriate, then this will be consistently estimated as well. It is not necessary for the correlation structure to be correct for the consistency of the marginal parameter estimates, at least when the correlation structure is fixed; indeed the common choice of independence is almost certain not to be appropriate. However the estimates of precision of the marginal parameter estimates do need to be adjusted to allow for the true correlation structure. This correction is done in the so-called “sandwich” estimator provided by the procedure.
The procedures have been written so that it is possible to fit models other than the standard ones. An important example of such a model is the application of the GEE methodology to ordinal categorical data. This application requires the data to be arranged in a particular form (as cummulative logits) and a particular correlation matrix (specified in _GEECORRELATION). The type of analyses are explained in Kenward et al. (1994) and the methodology described in that paper has been duplicated. Further details are given in Kenward & Smith (1995b).
An option (SCALEFACTOR) has been included that allows the user to decide whether or not the scale factor is fixed at its independence distributional default, or is estimated from the scaled residuals as in Liang & Zeger (1986), or is treated as a vector varying over time.
GEE calls three subsidiary procedures, _GEECODI, _GEECALLIN and _GEECALDIS to assist with the analysis. There is no need for the user to modify these procedures:
_GEECODI |
prints out the results of the iterative processes; |
|---|---|
_GEECALLIN |
calculates the fitted values and derivatives for various links; |
_GEECALDIS |
calculates the variance function and deviance for various distributions. |
There are also four other procedures, which can be re-written or replaced, to cater for further user-defined distributions, links and correlation structures:
_GEEINIT |
calculates initial estimates of the linear predictor in the generalized linear model; |
|---|---|
_GEELINK |
calculates fitted values and derivatives; |
_GEEDISTRIBUTION |
calculates the variance function and deviance; |
_GEECORRELATION |
calculates the correlation matrix and the sandwich matrix involving the residuals. (For the normal distribution the variance-covariance matrices are used not the correlation matrices.) |
If the LINK option is unset, the procedure will call _GEEINIT and _GEELINK instead of using those for the various standard link functions. For a logit link function _GEEINIT and _GEELINK should be defined as follows.
PROCEDURE ‘_GEEINIT’
“Calculation of initial estimate of linear predictor,
link unset”
PARAMETER NAME = \
‘Y’, “I: variate; response variate”\
‘LINEARPREDICTOR’,”O: variate; linear predictor”\
‘OFFSET’, “I: variate; offset”\
‘NBINOMIAL’; “I: variate; denominator of binomial”\
SET=3(yes),no;TYPE=4(‘variate’); \
COMPATIBLE=*,3(!T(type,nvalues,restriction));\
PRESENT=yes,no,2(yes)
CALC LINEARPREDICTOR = LOG((Y+0.5)/(NBINOMIAL-Y+0.5)) – OFFSET
ENDPROCEDURE
PROCEDURE ‘_GEELINK’
“Calculation of fitted values and derivatives”
PARAMETER NAME = \
‘LINEARPREDICTOR’, “I: variate; linear predictor”\
‘FITTEDVALUES’, “O: variate; estimate of fitted values”\
‘DERIVATIVES’, “O: variate; estimate of derivatives”\
‘OFFSET’, “I: variate; offset”\
‘NBINOMIAL’; “I: variate; denominator of binomial”\
SET=4(yes),no;TYPE=5(‘variate’); \
COMPATIBLE=*,4(!T(type,nvalues,restriction));\
PRESENT=yes,2(no),2(yes)
GETATTRIBUTE [ATTRIBUTE=NVALUES] LINEARPREDICTOR; SAVE=!P(nobs)
CALC FITTEDVALUES = NBINOMIAL/(1+EXP(-LINEARPREDICTOR – OFFSET))
& DERIVATIVES = 1/FITTEDVALUES+1/(NBINOMIAL-FITTEDVALUES)
ENDPROCEDURE
If the DISTRIBUTION option is unset, the procedure will call _GEEDISTRIBUTION instead of using one of the various standard distributions. For a binomial error distribution _GEEDISTRIBUTION should be defined as follows.
PROCEDURE ‘_GEEDISTRIBUTION’
” Calculation of variance function and deviance”
PARAMETER NAME = \
‘Y’, “I: variate; response variate”\
‘FITTEDVALUES’,”I: variate; fitted values”\
‘VARIANCE’, “O: variate; variance”\
‘DEVIANCE’, “O: scalar; total deviance”\
‘NBINOMIAL’; “I: variate; denominator of binomial”\
SET=4(yes),no;TYPE=3(‘variate’),’scalar’,’variate’; \
COMPATIBLE=*,2(!T(type,nvalues,restriction)),*,\
!T(type,nvalues,restriction); \
PRESENT=2(yes),2(no),yes
CALC VARIANCE = FITTEDVALUES*(NBINOMIAL-FITTEDVALUES)/NBINOMIAL
& DEVIANCE = -2*LLB(Y;NBINOMIAL;(FITTEDVALUES/NBINOMIAL))
ENDPROCEDURE
If the CRTYPE option is unset, the procedure will call _GEECORRELATION instead of using one of the various standard correlation models. For the independence model _GEECORRELATION should be defined as follows. Kenward & Smith (1995b) describe how _GEECORRELATION should be set up for analysing repeated ordinal categorical data.
PROCEDURE ‘_GEECORRELATION’
“
Calculation of correlation matrix
For SANDWICH = NO
input is the R matrix as for UNSPECIFIED
output is the desired R matrix.
For SANDWICH = YES
input is the (Y-MU)*T(Y-MU) matrix
output is the desired modified (Y-MU)*T(Y-MU) matrix.
N.B. For the normal distribution both the input and
output R’s should be variance/covariance matrices
not correlation matrices.
“
OPTION NAME = \
‘CONSTANT’, “I: text; how to treat constant (estimate,
omit); default e”\
‘SANDWICH’; “I; text; whether the sandwich central matrix
product or not) (no,yes); default no”\
MODE=2(T); NVALUES=2(1); SET=yes;\
VALUES=!T(ESTIMATE,OMIT),!T(NO,YES); \
DEFAULT=!T(ESTIMATE),!T(NO);
PARAMETER NAME = \
‘CORRELATIONS’,”I/O: matrix; the correlation matrix”\
ESTIMATES’, “I: variate; estimates of parameters in
model”\
‘Y’, “I: variate; response variate”\
‘RESIDUALS’, “I: variate; residuals”\
‘FITTEDVALUES’,”I: variate; fitted values”\
TIME’, “I: variate; times of repeated measures”\
‘MARKER’, “I: factor; identifier of subject or outcome”\
‘DISTRIBUTION’,”I: text; identifier of distribution”\
‘SCALEFACTOR’, “I: text; scalefactor option in use”\
‘SFVALUE’; “I: scalar; value of scalefactor if FIXED”\
SET=10(yes);DECLARED=10(yes); \
TYPE=’symmetric’,5(‘variate’),’factor’,2(‘text’),’scalar’; \
PRESENT=9(yes),no
GETATTRIBUTE [ATTRIBUTE=NVALUES] ESTIMATES; SAVE=!P(ncol)
& [ATTRIBUTE=NROWS] CORRELATIONS; SAVE=!P(ntime)
DIAGONALMATRIX [ROWS=ntime;MODIFY=yes] done,wkdm; \
VALUES=!(#ntime(1)),*
CALC const = ‘ESTIMATE’ .IN. CONSTANT
& sandw = ‘NO’ .IN. SANDWICH
IF sandw
“
SCALEFACTOR is as in GEE i.e. FIXED means fixed to SFVALUE
CONSTANT means the scalefactor is estimated but constant
across time, and VARYTIME means the scalefactor is estimated
and varies across time.
The variate TIME in this PROCEDURE represents the 1…ntime
distinct times, it is not a FACTOR of length nobs as in GEE.
It is the levels of the parameter TIME of GEE.
“
IF DISTRIBUTION.EQS.’NORMAL’
IF SCALEFACTOR.NES.’VARYTIME’
IF SCALEFACTOR.EQS.’FIXED’
CALC wkdm = SFVALUE
ELSE
CALC wkdm = TRACE(CORRELATIONS)/ntime
ENDIF
ELSE
CALC wkdm = CORRELATIONS
ENDIF
CALC CORRELATIONS = 0 + wkdm
ELSE
CALC CORRELATIONS = done
ENDIF
ENDIF
ENDPROCEDURE
If LINK, DISTRIBUTION or CRTYPE are unset, but no user routines are given for _GEEINIT, _GEELINK, _GEEDISTRIBUTION and _GEECORRELATION, then those given here (for logit link, binomial error distribution and independence) will be used.
Action with RESTRICT
Input structures must not be restricted, and any existing restrictions will be cancelled.
References
Kenward, M.G., Lesaffre, E. & Molenberghs, G. (1994). An application of maximum likelihood and generalized estimating equations to the analysis of ordinal data from a longitudinal study with cases missing at random. Biometrics, 50, 945-953.
Kenward, M.G. & Smith, D.M. (1995a). Computing the generalized estimating equations for repeated measurements. Genstat Newsletter, 32, 50-62.
Kenward, M.G. & Smith, D.M. (1995b). Computing the generalized estimating equations for repeated ordinal, categorical measurements. Genstat Newsletter, 32, 63-70.
Liang, K.-Y. & Zeger,S.L. (1986). Longitudinal data analysis using generalized linear models. Biometrika, 73, 13-22.
Liang, K.-Y., Zeger, S.L. & Qaqish, B. (1992). Multivariate regression analyses for categorical data (with discussion). Journal of the Royal Statistical Society, Series B, 54, 3-40.
See also
Commands for: Regression analysis, Repeated measurements.
Example
CAPTION 'GEE example',!t('Data from Archer',\
'(1987, Fertility & Sterility, 47, 559-564)',\
'about prolactin response to thyroptrin releasing',\
'hormone in women grouped according fertility status;',\
'also see Paik (1992, Biometrics, 48, 19-30).'); STYLE=meta,plain
VARIATE [VALUES=44,45,41,40,72,49,41,31,37,23,15,10,103,79,47,\
39,51,40,32,15,97,98,76,51,59,55,49,36,97,75,\
49,38,88,78,61,43,53,40,29,23,66,35,18,16,60,\
48,32,29,53,47,29,38,111,77,59,58,27,22,18,12,\
51,62,40,37,28,33,20,15,49,39,32,23,59,49,43,\
38,155,126,72,48,82,67,54,44,127,99,58,53,75,59,\
46,29,71,62,49,44,114,110,95,52,172,95,51,43,210,\
156,117,91,100,90,60,50,86,65,57,42,101,93,68,47] Response
& [VALUES=16,16,16,16,10,10,10,10,7,7,7,7,8,8,8,8,5,5,5,5,\
8,8,8,8,7,7,7,7,9,9,9,9,13,13,13,13,3,3,3,3,\
7,7,7,7,8,8,8,8,17,17,17,17,38,38,38,38,11,11,11,11,\
12,12,12,12,7,7,7,7,7,7,7,7,26,26,26,26,9,9,9,9,\
19,19,19,19,12,12,12,12,20,20,20,20,41,41,41,41,4,4,4,4,\
30,30,30,30,15,15,15,15,36,36,36,36,15,15,15,15,11,11,11,11]\
Baseline
& [VALUES=(1...4)30] CTime
FACTOR [LEVELS=30; VALUES=4(1...30)] Woman
& [LEVELS=3; VALUES=24(1),48(2),48(3)] Group
& [LEVELS=4; VALUES=(1...4)30] Time
GEE [LINK=log; DISTRIBUTION=gamma; TERMS=Group+CTime+Baseline;\
CRTYPE=AUTOREGRESSIVE] SUBJECT=Woman; TIME=Time; Y=Response