JSS

Journal of Statistical Software
March 2015, Volume 64, Issue 8.

doi: 10.18637/jss.v064.i08

R Package multgee: A Generalized Estimating
Equations Solver for Multinomial Responses
Anestis Touloumis
University of Brighton

Abstract
This vignette of the R package multgee is an updated version of Touloumis (2015),
published in the Journal of Statistical Software. To cite multgee in publications, please
use Touloumis (2015). To cite the GEE methodology implemented in multgee, please use
Touloumis, Agresti, and Kateri (2013).
The R package multgee implements the local odds ratios generalized estimating equations (GEE) approach proposed by Touloumis et al. (2013), a GEE approach for correlated
multinomial responses that circumvents theoretical and practical limitations of the GEE
method. A main strength of multgee is that it provides GEE routines for both ordinal(ordLORgee) and nominal (nomLORgee) responses, while relevant software in R and SAS
are restricted to ordinal responses under a marginal cumulative link model specification.
In addition, multgee offers a marginal adjacent categories logit model for ordinal responses
and a marginal baseline category logit model for nominal. Further, utility functions areavailable to ease the local odds ratios structure selection (intrinsic.pars), to constract
confidence intervals (confint), to extract the variance-covariance matrix of the regression
parameters (vcov) and to perform a Wald–type goodness-of-fit test between two nested
GEE models (waldts). We demonstrate the application of multgee through a clinical
trial with clustered ordinal multinomial responses.

Keywords: generalized estimating equations, nominal and ordinal multinomial responses, local
odds ratios, R.

1. Introduction
In several studies, the interest lies in drawing inference about the regression parameters ofa marginal model for correlated, repeated or clustered multinomial variables with ordinal or
nominal response categories while the association structure between the dependent responses
is of secondary importance. The lack of a convenient multivariate distribution for multinomial

2

multgee: GEE for Multinomial Responses

responses and the sensitivity of ordinary maximum likelihood methods to misspecification of
the association structure led researchers to modify the GEE method of Liang and Zeger (1986)
in order to account for multinomial responses (Miller, Davis, and Landis 1993; Lipsitz, Kim,
and Zhao 1994; Williamson, Kim, and Lipsitz 1995; Lumley 1996; Heagerty and Zeger 1996;
Parsons, Edmondson, and Gilmour 2006). These GEE approaches estimate the marginal
regression parameter vector by solving the same set of estimating equations as in Liang andZeger (1986), but differ in the way they parametrize and/or estimate α, a parameter vector
that is usually defined to describe a “working” assumption about the association structure.
Touloumis et al. (2013) showed that the joint existence of the estimated marginal regression parameter vector and α̂ cannot be assured in existing approaches. This is because theparametric space of the proposed parameterizations of the association structure depends onthe marginal model specification even in the simple case of bivariate multinomial responses.
To address this issue, Touloumis et al. (2013) defined α as a “nuisance” parameter vector
that contains the marginalized local odds ratios structure, that is the local odds ratios as if
no covariates were recorded, and they employed the family of association models (Goodman1985) to develop parsimonious and meaningful structures regardless of the response scale.
The practical advantage of the local odds ratios GEE approach is that it is applicable to bothordinal and nominal multinomial responses without being restricted by the marginal model
specification. Simulations in Touloumis et al. (2013) imply that the local odds ratios GEE
approach captures a significant portion of the underlying correlation structure, and comparedto the independence “working” model (i.e., assuming no correlation structure in the GEE
methodology), simple local odds ratios structures can substantially increase the efficiency
gains in estimating the regression vector of the marginal model. Note that low convergence
rates for the GEE approach of Lipsitz et al. (1994) and Heagerty and Zeger (1996) did notallow the authors to compare these approaches with the local odds ratios GEE approach while
the GEE approach of Parsons et al. (2006) was excluded from the simulation design because
its use is restricted to a cumulative logit marginal model specification.
The R (R Core Team 2019) package multgee implements the local odds ratios GEE approachand it is available from the Comprehensive R Archive Network at https://CRAN.R-project.
org/package=multgee. To emphasize the importance of reflecting the nature of the response
scale on the marginal model specification and on the marginalized local odds ratios structure,
two core functions are available in multgee: nomLORgee which is appropriate for GEE analysis
of nominal multinomial responses and ordLORgee which is appropriate for ordinal multinomial
responses. In particular, options for the marginal model specification include a baseline
category logit model for nominal response categories and a cumulative link model or an
adjacent categories logit model for ordinal response categories. In addition, there are six
utility functions that enable the user to: i) Perform goodness-of-fit tests between two nested
GEE models (waldts), ii) construct Wald type confidence intervals (confint), iii) select the
local odds ratios structure based on the rule of thumb discussed in Touloumis et al. (2013)(intrinsic.pars), iv) obtain the estimated covariance matrix of the parameter estimates(vcov), v) construct a probability table (to be passed in the core functions) that satisfies adesired local odds ratios structure (matrixLOR), and vi) obtain commonly used criteria forvariables and/or association structure selection in the GEE literature.
To appreciate the features of multgee, we briefly review GEE software for multinomial responses in SAS (SAS Institute Inc. 2014) and R. The current version of SAS supports only the
independence “working” model under a marginal cumulative probit or logit model for ordinal

Journal of Statistical Software3

multinomial responses. To the best of our knowledge, SAS macros (Williamson, Lipsitz, andKim 1998; Yu and Yuan 2004) implementing the approach of Williamson et al. (1995) are
not publicly available. The R package repolr (Parsons 2013) implements the approach of
Parsons et al. (2006) but it is restricted to using a cumulative logit model. Another option
for ordinal responses is the function ordgee in the R package geepack (Halekoh, Højsgaard,
and Yan 2006). This function implements the GEE approach of Heagerty and Zeger (1996)
but it seems to produce unreliable results for multinomial responses. To illustrate this, we
simulated independent multinomial responses under a cumulative probit model specification
with a single time-stationary covariate for each subject and we employed ordgee to obtain
the GEE estimates from the independence “working” model. Description of the generative
process can be found in Scenario 1 of Touloumis et al. (2013) except that we used the values−3, −1, 1 and 3 for the four category specific intercepts in order to make the problem moreevident. Based on 1000 simulation runs when the sample size N = 500, we found that thebias of the GEE estimate of β = 1 was ≈ 4.8 × 1028 , indicating the presence of a bug or -at
least- of numerical problems for some situations. Similar problems occurred for the alternative global odds ratios structures in ordgee. In contrast to existing software, multgee offers
greater variety of GEE models for ordinal responses, implements a GEE model for nominal
responses and is not limited to the independence “working” model, which might lead to significant efficiency losses. Further, one can assess the goodness of fit for two or more nested
GEE models.
This paper is organized as follows. In Section 2, we present the theoretical background of
the local odds ratios GEE approach that is necessary for the use of multgee. We introduce
the marginal models implemented in multgee, the estimation procedure for the “nuisance”
parameter vector α and the asymptotic theory on which GEE inference is based. We describe
the arguments of the core GEE functions (nomLORgee, ordLORgee) in Section 3 while theutility functions (waldts, confint, intrinsic.pars, vcov, matrixLOR, gee_criteria) aredescribed in Section 4. In Section 5, we illustrate the use of multgee in a longitudinal study
with correlated ordinal multinomial responses. We summarize the features of the package
and provide a few practical guidelines in Section 6.

2. Local odds ratios GEE approach
For notational ease, suppose the data arise from a longitudinal study with no missing observations. However, note that the local odds ratios GEE approach is not limited neither
to longitudinal studies nor to balanced designs, under the strong assumption that missing
observations are missing completely at random (Rubin 1976).
Let Yit be the multinomial response for subject i (i = 1, . . . , N ) at time t (t = 1, . . . , T ) thattakes values in ¶1, 2, . . . , J♢, J > 2. Define the response vector for subject i
Yi = (Yit1 , . . . , Yi1(J−1) , Yi21 , . . . , Yi2(J−1) , . . . , YiT 1 , . . . , YiT (J−1) )⊤ ,where Yitj = 1 if the response for subject i at time t falls at category j and Yitj = 0 otherwise.
⊤ , . . . , x⊤ )⊤ be theDenote by xit the covariates vector associated with Yit , and let xi = (xi1iT
covariates matrix for subject i. Define πitj = E(Yitj ♣xi ) = P(Yitj = 1♣xi ) = P(Yit = j♣xi ) as
⊤ ⊤the probability of the response category j for subject i time t, and let π i = (π ⊤i1 , . . . , π iT )⊤be the mean vector of Yi , where π it = (πit1 , . . . , πit(J−1) ) . It follows from the above thatP
PJ−1
YitJ = 1 − J−1
j=1 Yitj and πitJ = 1 −j=1 πitj .

4multgee: GEE for Multinomial Responses

2.1. Marginal models for correlated multinomial responses
The choice of the marginal model depends on the nature of the response scale. For ordinal
multinomial responses, the family of cumulative link models
F −1 [P(Yit ≤ j♣xi )] = βj0 + β ⊤
∗ xit

(1)

or the adjacent categories logit model
log



πitj
πit(j+1)= βj0 + β ⊤
∗ xit

(2)

can be used, where F is the cumulative distribution function of a continuous distributionand ¶βj0 : j = 1, . . . , J − 1♢ are the category specific intercepts. For nominal multinomial
responses, the baseline category logit model
log



πitj
πitJ



= βj0 + β ⊤
j xit

(3)

can be used, where β j is the j-th category specific parameter vector.
It is worth mentioning that the linear predictor differs in the above marginal models. First, thecategory specific intercepts need to satisfy a monotonicity condition β10 ≤ β20 ≤ . . . ≤ β(J−1)0
only when the family of cumulative link models in (1) is employed. Second, the regression
parameter coefficients of the covariates xit are category specific only in the marginal baseline
category logit model (3) and not in the ordinal marginal models (1) and (2).

2.2. Estimation of the marginal regression parameter vector
To unify the notation, let β be the p–variate parameter vector that includes all the regression
parameters in (1), (2) or (3). To obtain βbG , a GEE estimator of β, Touloumis et al. (2013)
solved the estimating equationsb =U(β, α)N
1 X
Di Vi−1 (Yi − π i ) = 0
N i=1

(4)

where Di = ∂π i /∂β and Vi is a T (J − 1) × T (J − 1) “weight” matrix that depends on βb an estimate of the “nuisance” parameter vector α defined formally in Section 2.3.
and on α,
Succinctly, Vi is a block matrix that mimics the form of COV(Yi ♣xi ), the true covariancematrix for subject i. The t-th diagonal matrix of Vi is the covariance matrix of Yit determined by the marginal model. The (t, t′ )-th off-diagonal block matrix describes the marginal
pseudo-association of (Yit , Yit′ ), which is a function of the marginal model and of the pseudob β). We
probabilities ¶P(Yit = j, Yit′ = j ′ ♣xi ) : j, j ′ = 1, . . . , J − 1♢ calculated based on (α,should emphasize that Vi is a “weight” matrix because α is defined as a “nuisance” parameter vector and it is unlikely to describe a valid “working” assumption about the association
structure for all subjects.

2.3. Estimation of the nuisance parameter vector and of the weight matrix
Order the L = T (T − 1)/2 time-pairs with the rightmost element of the pair most rapidly
varying as (1, 2), (1, 3), . . . , (T − 1, T ), and let G be the group variable with levels the L

5Journal of Statistical Softwareordered pairs. For each time-pair (t, t′ ), ignore the covariates and cross-classify the responses
across subjects to form an J × J contingency table such that the row totals correspond tothe observed totals at time t and the column totals to the observed totals at time t′ , andlet θtjt′ j ′ be the local odds ratio at the cutpoint (j, j ′ ) based on the expected frequencies¶ftjt′ j ′ : j, j ′ = 1, . . . , J♢. For notational reasons, let A and B be the row and column variable
respectively. Assuming a Poisson sampling scheme to the L sets of J × J contingency tables,
fit the RC-G(1) type model (Becker and Clogg 1989)(t,t′ ) (t,t′ )µj ′ ,′(t,t )BG
AG
G
Bµj
log ftjt′ j ′ = λ + λA
j + λj ′ + λ(t,t′ ) + λj(t,t′ ) + λj ′ (t,t′ ) + ϕ(5)(t,t′ )

where ¶µj
: j = 1, . . . , J♢ are the score parameters for the J response categories at the′
time-pair (t, t ). After imposing identifiability constraints on the regression parameters in (5),the log local odds ratios structure is given by′(t,t′ )log θtjt′ j ′ = ϕ(t,t ) µj(t,t′ )

− µj+1(t,t′ )

µj ′(t,t′ )

− µj ′ +1 .

(6)

At each time-pair, (6) summarizes the local odds ratios structure in terms of the J score′parameters and the intrinsic parameter ϕ(t,t ) that measures the average association of themarginalized contingency table. Since the score parameters do not need to be fixed or monotonic, the local odds ratios structure is applicable to both nominal and ordinal multinomial
responses.
Touloumis et al. (2013) defined α as the parameter vector that contains the marginalized
local odds ratios structure

α = θ1121 , . . . , θ1(J−1)2(J−1) , . . . , θ(T −1)1T 1 , . . . , θ(T −1)(J−1)T (J−1)

⊤where θtjt′ j ′ satisfy (6). To increase the parsimony of the local odds ratios structuresfor
(t,t′ )
ordinal responses, they proposed to use common unit-spaced score parameters µj
=j


′



and/or common intrinsic parameters ϕ(t,t ) = ϕ across time-pairs. For a nominal response(t,t′ )

scale, they proposed to apply a homogeneity constraint on the score parameters µj= µj
and use common intrinsic parameters across time-pairs. To estimate α maximum likelihood
methods are involved by treating the L marginalized contingency tables as independent.
Technical details and justification about this estimation procedure can be found in Touloumis(2011) and Touloumis et al. (2013).
b and the marginalConditional on the estimated marginalized local odds ratios structure α
model specification at times t and t′ , ¶P(Yit = j, Yit′ = j ′ ♣xi ) : t < t′ , j, j ′ = 1, . . . , J − 1♢ are
obtained as the unique solution of the iterative proportional fitting (IPF) procedure (Deming
and Stephan 1940). Hence, Vi can be readily calculated and the estimating equations in (4)can be solved with respect to β.

2.4. Asymptotic properties of the GEE estimator√b − β) ∼ N (0, Σ) asymptotically,b inference about β is based on the fact that N (βGiven α,Gwhere−1
Σ = lim N Σ−1(7)
0 Σ1 Σ0 ,N →∞6multgee: GEE for Multinomial Responses
P

P

−1
N
⊤ −1
⊤ −1
Σ0 = Ni=1 Di Vi Di and Σ1 =i=1 Di Vi COV(Yi ♣xi )Vi Di . For finite sample sizes,b and COV(Yi ♣xi ) withb is estimated by ignoring the limit in (7) and replacing β with βΣG⊤bb i )(Yi − πb i ) in Σ0 and Σ1 . In the literature, Σ/N(Yi − πis often termed as “sandwich” or“robust” covariance matrix of βbG .

3. Description of core functions
We describe the arguments of the functions nomLORgee and ordLORgee, focusing on themarginal model specification (formula, link), data representation (id, repeated, data) andlocal odds ratios structure specification (LORstr, LORterm, homogeneous, restricted). For
completeness’ sake, we also present computational related arguments (LORem, add, bstart,
LORgee_control, ipfp.control, IM). The two core functions share the same arguments,
except link and restricted which are available only in ordLORgee, and they both create an object of the class LORgee which admits summary, coef, update, vcov, confint andresiduals methods.

3.1. Marginal model specification
For ordinal multinomial responses, the link argument in the function ordLORgee specifies
which of the marginal models (1) or (2) will be fitted. The options "logit", "probit","cauchit" or "cloglog" indicate the corresponding cumulative distribution function F in
the cumulative link model (1), while the option "acl" implies that the adjacent categories
logit model (2) is selected. For nominal multinomial responses, the function nomLORgee fits
the baseline category logit model (3), and hence the link argument is not offered.
The formula (= response ~ covariates) argument identifies the multinomial response variable (response) and specifies the form of the linear predictor (covariates), assuming that
this includes an intercept term. If required, the J > 2 observed response categories are sorted
in an ascending order and then mapped onto ¶1, 2, . . . , J♢. To account for a covariate x witha constrained parameter coefficient fixed to 1 in the linear predictor, the term offset(x)must be inserted on the right hand side of formula.

3.2. Data representation
The id argument identifies the N subjects by assigning a unique label to each subject. If required, the observed id labels are sorted in an ascending order and then relabeled as 1, . . . , N ,respectively.
The repeated argument identifies the times at which the multinomial responses are recordedby treating the T unique observed times in the same manner as in id. The purpose of
repeated is dual: To identify the T distinct time points and to construct the full marginalized
contingency table for each time-pair by aggregating the relevant/available responses acrosssubjects. The repeated argument is optional and it can be safely ignored in balanced designs
or in unbalanced designs in which if the t-th response is missing for a particular subject
then all subsequent responses at times t′ > t are missing for that subject. Otherwise, it isrecommended to provide the repeated argument in order to ensure proper construction of
the full marginalized contingency table. To this end, note that if the measurement occasions
are not recorded in a numerical mode, then the user should create repeated by mapping the7

Journal of Statistical Softwarelog θtjt′ j ′
ϕ
′
(t,t
ϕ )

ϕ
(µj − µj+1 ) µj′ − µj ′ +1
′(t,t′ )
(t,t′ )
(t,t′ )
(t,t′ )
µj ′ − µj ′ +1
− µj+1ϕ(t,t ) µj

LORstr"uniform""category.exch""time.exch""RC"

Functions
ordLORgee
ordLORgee
Both
Both

Parameters
1
L
J −1
L(J − 1)

Table 1: The main options for the marginalized local odds ratios structures in multgee.
T distinct measurement occasions onto the set ¶1, . . . , T ♢ in such a way that the temporalorder of the measurement occasions is preserved. For example, if the measurements occasions
are recorded as “before”, “baseline”, “after”, then the levels for repeated should be coded as
1, 2 and 3, respectively.
The dataset is imported via the data argument in “long” format, meaning that each row
contains all the information provided by a subject at a given measurement occasion. Thisimplies that data must include the variables specified in the mandatory arguments formulaand id, as well as the optional argument repeated when this is specified by the user. If nodata is provided then the above variables are extracted from the environment that nomLORgee
and ordLORgee are called. Currently missing observations, identified by NA in data, areignored.

3.3. Marginalized local odds ratios structure specification
The marginalized local odds ratios is specified via the LORstr argument. Table 1 displays
the structures proposed by Touloumis et al. (2013). Currently the default option is thetime exchangeability structure ("time.exch") in nomLORgee and the category exchangeability ("category.exch") structure in ordLORgee. The uniform ("uniform") and category exchangeability structures are not allowed in nomLORgee because given unit-spaced parameterscores are not meaningful for nominal response categories.
The user can also fit the independence “working” model (LORstr = "independence") or evenprovide the local odds ratios structure (LORstr = "fixed") using the LORterm argument. In
this case, an L×J 2 matrix must be constructed such that the g-th row contains the vectorized
form of a probability table that satisfies the desired local odds ratios structure at the time-paircorresponding to the g-th level of G.
Touloumis (2011) discussed two further versions of the "time.exch" and the RC ("RC")structures based on using: i) Heterogeneous score parameters (homogeneous = FALSE) at eachtime-pair, and/or ii) monotone score parameters (restricted = TRUE), an option applicable
only for ordinal response categories. However, it is sensible to employ these additional options
only when the local odds ratios structures in Table 1 do not seem adequate.
It is important to mention that the user must provide only the arguments required for the
specified local odds ratios structure. For example, the arguments homogeneous, restricted
and LORterm are ignored when LORstr = "uniform".

3.4. Computational details
The default estimation procedure for the marginalized local odds ratios structure is to fit
model (5) to the full marginalized contingency table (LORem = "3way") after imposing the8multgee: GEE for Multinomial Responses

desired restrictions on the intrinsic and the score parameters. Touloumis (2011) noticed thatthe estimated local odds ratios structure under model (5) is identical to that obtained by fitting independently a row and columns (RC) effect model (Goodman 1985) with homogeneousscore parameters to each of the L contingency tables. Motivated by this, an alternative estimation procedure (LORem = "2way") for estimating the structures "uniform" and "time.exch"
was proposed. In particular, one can estimate the single parameter of the "uniform" struc′ture as the average of the L intrinsic parameters ϕ(t,t ) obtained by fitting the linear-by-linear
association model (Agresti 2013) independently to each of the L marginalized contingency
tables. For the "time.exch" structure, one can fit L RC effects models with homogeneous(homogeneous = TRUE)/heterogeneous (homogeneous = FALSE) score parameters and thenestimate the log local odds ratio at each cutpoint (j, j ′ ) by averaging log θ̂tjt′ j ′ for t < t′ .
Regardless of the value of LORem, the appropriate model for counts is fitted via the function
gnm of the R package gnm (Turner and Firth 2012).
In the presence of zero observed counts, a small positive constant can be added (add) at eachb We conjecture that acell of the marginalized contingency table to ensure the existence of α.
constant of the magnitude 10−4 will serve this purpose without affecting the strength of theassociation structure.
A Fisher scoring algorithm is employed to solve the estimating equations (4) as in Lipsitz et al.
b is not updated. The default way to obtain the initial(1994). The only difference is that now α
value for β is via the function vglm of the R package VGAM (Yee 2010). Alternatively, theinitial value can be provided by the user (bstart). The Fisher scoring algorithm converges
when the elementwise maximum relative change in two consecutive estimates of β is lessthan or equal to a predefined positive constant ϵ. The control argument controls the related
iterative procedure variables and printing options. The default maximum number of iterations
is 15 and the default tolerance is ϵ = 0.001.
Recall that calculation of the “weight” matrix Vi at given values of (β, α) relies on the IPF
procedure. The ipfp.ctrl argument controls the related variables. The convergence criterion
is the maximum of the absolute difference between the fitted and the target row and column
marginals. By default, the tolerance of the IPF procedure is 10−6 with a maximal number of
iterations equal to 200.
The IM argument defines which of the R functions solve, qr.solve or cholesky will be used
to invert matrices in the Fisher scoring algorithm.

4. Description of utility functions
The function waldts performs a goodness-of-fit test for two nested GEE models based on a
Wald test statistic. Let M0 and M1 be two nested GEE models with marginal regression pa⊤ ⊤rameter vectors β 0 and β 1 = (β ⊤0 , β q ) , respectively. Define a matrix C such that Cβ 1 = β q .
Here q equals the rank of C and the dimension of β q . The hypothesis
H0 : β q = 0 vs H1 : β q ̸= 0
tests the goodness-of-fit of M0 versus M1 . Based on a Wald type approach, H0 is rejectedb ⊤ (N CΣCb ≥ Xq (α), where βb and Σb ⊤ )−1 (Cβ)
b are estimated
at α% significance level, if (Cβ)
under model M1 and Xq (α) denotes the α upper quantile of a chi-square distribution with q
degrees of freedom.

Journal of Statistical Software

9

The function confint calculates the Wald-type confidence intervals for the regression parameters of the marginal model under consideration. The method argument specifies whether the
standard errors of the regression parameters estimates will be derived from the “sandwich”(method = "robust") or the “model-based” (method = "naive") covariance matrix. The
default is to use the “sandwich” covariance matrix.
Touloumis et al. (2013) suggested to select the local odds ratios structure by inspecting
the range of the L estimated intrinsic parameters under the "category.exch" structure forordinal responses, or under the "RC" structure for nominal responses. If the estimated intrinsic
parameters do not differ much, then the underlying marginalized local odds ratios structure is
likely nearly exchangeable across time-pairs. In this case, the simple structures "uniform" or"time.exch" should be preferred because they tend to be as efficient as the more complicatedones. The function intrinsic.pars gives the estimated intrinsic parameter of each time-pair.
The function vcov returns the “sandwich” (method = "robust") or the “model-based” (method= "naive") covariance matrix. The default is to use the “sandwich” covariance matrix. It is
used internally to calculate confidence intervals for the regression parameters and to compare
nested models.
The single-argument function matrixLOR creates a two-way probability table that satisfies adesired local odds ratios structure. This function aims to ease the construction of the LORterm
argument in the core functions nomLORgee and ordLORgee.
For completeness, the function gee_criteria reports commonly used criteria for variableselection and for selecting the “working” association structure. The Quasi Information Criterion (QIC) proposed by Pan (2001), the Correlation Information Criterion (CIC) proposed by
Hin and Wang (2009) and the Rotnitzky and Jewell’s Criterion (RJC) proposed by Rotnitzky
and Jewell (1990) are used for selecting the best association structure. The QICu criterion(Pan 2001) is used for selecting the best subset of covariates. When choosing among GEE
models with different association structures but with the same subset of covariates, the modelwith the smallest value of QIC, CIC or RJC should be preferred. However, note that thesecriteria might not be very useful for multinomial GEE models with local odds ratios structures and the user should prefer the empirical rule suggested in Touloumis et al. (2013) and
implemented in intrinsic.pars. When choosing between two or more GEE models withdifferent number of covariates that are not nested, the model with the smallest QICu value
should be preferred. The function gee_criteria can be used to obtain the values for theabove criteria for one or more models of the class LORgee.

5. Example
To illustrate the main features of the package multgee, we follow the GEE analysis performedin Touloumis et al. (2013). The data came from a randomized clinical trial (Lipsitz et al. 1994)
that aimed to evaluate the effectiveness of the drug Auranofin versus the placebo therapy forthe treatment of rheumatoid arthritis. The five-level (1=poor, . . . , 5=very good) ordinal
multinomial response variable was the self-assessment of rheumatoid arthritis recorded at one(t = 1), three (t = 2) and five (t = 3) follow-up months. To acknowledge the ordinal response

10

multgee: GEE for Multinomial Responses

scale, the marginal cumulative logit modelP(Yit ≤ j♣xi )log
1 − P(Yit ≤ j♣xi )

= βj0 + β1 I(timei = 3) + β2 I(timei = 5) + β3 trti
+ β4 I(bi = 2) + β5 I(bi = 3) + β6 I(bi = 4) + β7 I(bi = 5).

(8)

was fitted, where i = 1, . . . , 301, t = 1, 2, 3, j = 1, 2, 3, 4 and I(A) is the indicator function
for the event A. Here xi denotes the covariates matrix for subject i that includes the selfassessment of rheumatoid arthritis at the baseline (bi ), the treatment variable (trti ), codedas (1) for the placebo group and (2) for the drug group, and the follow-up time recorded inmonths (timei ).
The GEE analysis is performed in two steps. First, we select the marginalized local odds
ratios structure by estimating the intrinsic parameters under the "category.exch" structure
R> library("multgee")
Loading required package: gnm
R> data("arthritis")
R> head(arthritis)

1
2
3
4
5
6

id y sex age trt baseline time
1 4
2 54
2
2
1
1 5
2 54
2
2
3
1 5
2 54
2
2
5
2 4
1 41
1
3
1
2 4
1 41
1
3
3
2 4
1 41
1
3
5

R> intrinsic.pars(y = y, data = arthritis, id = id, repeated = time,+
rscale = "ordinal")[1] 0.6517843 0.9097341 0.9022272
The range of the estimated intrinsic parameters is small (≈ 0.26) which suggests that theunderlying marginalized association pattern is nearly constant across time-pairs. Thus we
expect the "uniform" structure to capture adequately the underlying correlation pattern.
Note that we passed the time variable to the repeated argument because this numericalvariable indicates the measurement occasion at which each observation was recorded.
Now we fit the cumulative logit model (8) under the "uniform" via the function ordLORgee
R>
+
+
R>

fit <- ordLORgee(formula = y ~ factor(time) + factor(trt) + factor(baseline),link = "logit", id = id, repeated = time, data = arthritis,
LORstr = "uniform")summary(fit)11

Journal of Statistical Software
GEE FOR ORDINAL MULTINOMIAL RESPONSES
version 1.6.0 modified 2017-07-10
Link : Cumulative logit
Local Odds Ratios:
Structure:
uniform
Model:
3way
call:
ordLORgee(formula = y ~ factor(time) + factor(trt) + factor(baseline),
data = arthritis, id = id, repeated = time, link = "logit",
LORstr = "uniform")
Summary of residuals:Min.
1st Qu.
Median-0.5161111 -0.2399227 -0.0749670

Mean
3rd Qu.
0.0000219 -0.0066995

Max.
0.9932912

Number of Iterations: 5
Coefficients:
Estimatesan.sesan.z Pr(>|san.z|)
beta10-1.84315 0.38929 -4.7346< 2e-16 ***beta200.26692 0.35013 0.7624
0.44585
beta302.23132 0.36625 6.0924< 2e-16 ***beta404.52542 0.42123 10.7434< 2e-16 ***factor(time)3
0.00140 0.12183 0.0115
0.99080
factor(time)5-0.36172 0.11395 -3.1743
0.00150 **factor(trt)2-0.51212 0.16799 -3.0486
0.00230 **factor(baseline)2 -0.66963 0.38036 -1.7605
0.07832 .
factor(baseline)3 -1.26070 0.35252 -3.5763
0.00035 ***factor(baseline)4 -2.64373 0.41282 -6.4041< 2e-16 ***factor(baseline)5 -3.96613 0.53164 -7.4602< 2e-16 ***--Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1
Local Odds Ratios Estimates:[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12][1,] 0.000 0.000 0.000 0.000 2.257 2.257 2.257 2.257 2.257 2.257 2.257 2.257[2,] 0.000 0.000 0.000 0.000 2.257 2.257 2.257 2.257 2.257 2.257 2.257 2.257[3,] 0.000 0.000 0.000 0.000 2.257 2.257 2.257 2.257 2.257 2.257 2.257 2.257[4,] 0.000 0.000 0.000 0.000 2.257 2.257 2.257 2.257 2.257 2.257 2.257 2.257[5,] 2.257 2.257 2.257 2.257 0.000 0.000 0.000 0.000 2.257 2.257 2.257 2.257[6,] 2.257 2.257 2.257 2.257 0.000 0.000 0.000 0.000 2.257 2.257 2.257 2.257[7,] 2.257 2.257 2.257 2.257 0.000 0.000 0.000 0.000 2.257 2.257 2.257 2.257[8,] 2.257 2.257 2.257 2.257 0.000 0.000 0.000 0.000 2.257 2.257 2.257 2.257[9,] 2.257 2.257 2.257 2.257 2.257 2.257 2.257 2.257 0.000 0.000 0.000 0.000

12multgee: GEE for Multinomial Responses[10,] 2.257 2.257 2.257 2.257 2.257 2.257 2.257 2.257 0.000 0.000 0.000 0.000[11,] 2.257 2.257 2.257 2.257 2.257 2.257 2.257 2.257 0.000 0.000 0.000 0.000[12,] 2.257 2.257 2.257 2.257 2.257 2.257 2.257 2.257 0.000 0.000 0.000 0.000
p-value of Null model: < 0.0001
The summary method summarizes the fit of the GEE model including the GEE estimates,
their estimated standard errors based on the “sandwich” covariance matrix and the p–values
from testing the statistical significance of each regression parameter in (8). The estimatedmarginalized local odds ratios structure can be found in a symmetric T (J − 1) × T (J − 1)block matrix written symbolically as


0

 Θ21


 ...

Θ12
0

...
ΘT 1 ΘT 2

. . . Θ1T
. . . Θ2T .
..
. ... ...
0Each block denotes an (J −1)×(J −1) matrix. The (j, j ′ )-th element of the off-diagonal blockΘtt′ represents the estimate of θtjt′ j ′ . Based on the properties of the local odds ratios it is easy′
to see that Θtt′ = Θ⊤t′ t for t < t . Finally, the diagonal blocks are zero to reflect the fact thatno local odds ratios are estimated when t = t′ . In our example, J = 5 and thus each block
is a 4 × 4 matrix. Since the uniform structure is selected, all local odds ratios are equal andestimated as 2.257. Finally, pvalue of Null model corresponds to the p–value of testing
the hypothesis that no covariate is significant, i.e., β1 = β2 = β3 = β4 = β5 = β6 = β7 = 0,
based on a Wald test statistic.
The goodness-of-fit of model (8) can be tested by comparing it to a marginal cumulative logit
model that additionally contains the age and gender main effects in the linear predictor
R> fit1 <- update(fit, formula = ~. + factor(sex) + age)
R> waldts(fit, fit1)
Goodness of Fit based on the Wald test
Model under H_0: y ~ factor(time) + factor(trt) + factor(baseline)
Model under H_1: y ~ factor(time) + factor(trt) + factor(baseline) + factor(sex) +age
Wald Statistic = 3.9554, df = 2, p-value = 0.1384
Alternatively, we can compare these two models using commonly used criteria:
R> gee_criteria(fit, fit1)
QIC
CIC
RJC
QICu Parameters
fit 1264.325 1264.312 11.0066 1.2740
11
fit1 1275.860 1275.703 13.0785 0.9785
13

Journal of Statistical Software

13

The Wald confidence intervals for the regression parameters of the reduced model are:
R> confint(fit)
2.5 %
97.5 %
beta10-2.6061497 -1.08015721
beta20-0.4193163 0.95314883
beta301.5134867 2.94915281
beta403.6998267 5.35101277
factor(time)3-0.2373805 0.24018865
factor(time)5-0.5850535 -0.13837735
factor(trt)2-0.8413778 -0.18287050
factor(baseline)2 -1.4151270 0.07587068
factor(baseline)3 -1.9516257 -0.56978096
factor(baseline)4 -3.4528305 -1.83462340
factor(baseline)5 -5.0081187 -2.92413445

6. Summary and practical guidelines
We described the R package multgee which implements the local odds ratios GEE approach(Touloumis et al. 2013) for correlated multinomial responses. Unlike existing GEE softwares,
multgee allows GEE models for ordinal (ordLORgee) and nominal (nomLORgee) responses.
The available local odds ratios structures (LORstr) in each function respect the nature of the
response scale to prevent usage of ordinal local odds ratios structures (e.g., "uniform") innomLORgee. The fitted GEE model is summarized via the summary method while the estimatedregression coefficient can be retrieved via the coef method. The statistical significance of theregression parameters can be assessed via the function waldts. Confidence intervals of adesired marginal model can be obtained by the function confint. A similar strategy tothat presented in Section 5, can be adopted to analyze GEE models for correlated nominal
multinomial responses.
From a practical point of view, we recommend the use of the "uniform" structure for ordinal responses and the "time.exch" structure for nominal especially when the range of theestimated intrinsic parameters (intrinsic.pars) is small. Based on our experience, someconvergence problems might occur as the complexity of the local odds ratios structure increases and/or if the marginalized contingency tables are very sparse. Two possible solutions
are either to adopt a simpler local odds ratios structure or to increase slightly the value ofthe constant added to the marginalized contingency tables (add). However, we believe thatusers should refrain from using the independence “working” model unless the aforementioned
strategies fail to remedy the convergence problems. To decide on the form of the linear predictor, variable selection model procedures could be incorporated using the function waldts.
Finally, confidence intervals could be obtained by using the function confint.
In future versions of multgee, we plan to permit time-dependent intercepts in the marginalmodels, to increase the range of the marginal models, by including, for example, the familyof continuation-ratio models for ordinal responses, and to offer a function for assessing the
proportional odds assumption in models (1) and (2).

14

multgee: GEE for Multinomial Responses

References
Agresti A (2013). Categorical Data Analysis. 3rd edition. John Wiley & Sons.
Becker M, Clogg C (1989). “Analysis of Sets of Two-Way Contingency Tables Using Association Models.” Journal of the American Statistical Association, 84, 142–151.
Deming W, Stephan F (1940). “On a Least Squares Adjustment of a Sampled Frequency Table
When the Expected Marginal Totals Are Known.” The Annals of Mathematical Statistics,
11, 427–444.
Goodman L (1985). “The Analysis of Cross-Classified Data Having Ordered and/or Unordered Categories: Association Models, Correlation Models, and Asymmetry Models forContingency Tables With or Without Missing Entries.” The Annals of Statistics, 13, 10–69.
Halekoh U, Højsgaard S, Yan J (2006). “The R Package geepack for Generalized Estimating
Equations.” Journal of Statistical Software, 15, 1–11.
Heagerty P, Zeger S (1996). “Marginal Regression Models for Clustered Ordinal Measurements.” Journal of the American Statistical Association, 91, 1024–1036.
Hin LY, Wang YG (2009). “Working-correlation-structure identification in generalized estimating equations.” Statistics in medicine, 28(4), 642–658.
Liang K, Zeger S (1986). “Longitudinal Data Analysis Using Generalized Linear Models.”
Biometrika, 73, 13–22.
Lipsitz S, Kim K, Zhao L (1994). “Analysis of Repeated Categorical Data Using Generalized
Estimating Equations.” Statistics in Medicine, 13, 1149–1163.
Lumley T (1996). “Generalized Estimating Equations for Ordinal Data: A Note on the
Working Correlation Structures.” Biometrics, 52, 354–361.
Miller M, Davis C, Landis J (1993). “The Analysis of Longitudinal Polytomous Data: Generalized Estimating Equations and Connections with Weighted Least Squares.” Biometrics,
49, 1033–1044.
Pan W (2001). “Akaike’s information criterion in generalized estimating equations.” Biometrics, 57(1), 120–125.
Parsons N (2013). repolr: Repeated Measures Proportional Odds Logistic Regression. R
package version 2.0, URL https://CRAN.R-project.org/package=repolr.
Parsons N, Edmondson R, Gilmour S (2006). “A Generalized Estimating Equation Method for
Fitting Autocorrelated Ordinal Score Data with an Application in Horticultural Research.”
Journal of the Royal Statistical Society C, 55, 507–524.
R Core Team (2019). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/.
Rotnitzky A, Jewell NP (1990). “Hypothesis testing of regression parameters in semiparametric generalized linear models for cluster correlated data.” Biometrika, 77(3), 485–497.

15

Journal of Statistical Software
Rubin D (1976). “Inference and Missing Data.” Biometrika, 63, 581–592.

SAS Institute Inc (2014). SAS/STAT Software, Version 9.4. Cary, NC. URL http://www.
sas.com/.
Touloumis A (2011). Generalized Estimating Equations for Multinomial Responses. Ph.D.
thesis, University of Florida.
Touloumis A (2015). “R Package multgee: A Generalized Estimating Equations Solver for
Multinomial Responses.” Journal of Statistical Software, 64, 1–14.
Touloumis A, Agresti A, Kateri M (2013). “Generalized Estimating Equations for Multinomial
Responses Using a Local Odds Ratio Parameterization.” Biometrics, 69, 633–640.
Turner H, Firth D (2012). Generalized Nonlinear Models in R: An Overview of the gnm
Package. R package version 1.0-7, URL https://CRAN.R-project.org/package=gnm.
Williamson J, Kim K, Lipsitz S (1995). “Analyzing Bivariate Ordinal Data Using a Global
Odds Ratio.” Journal of the American Statistical Association, 90, 1432–1437.
Williamson J, Lipsitz S, Kim K (1998). “GEECAT and GEEGOR: Computer Programs for the
Analysis of Correlated Categorical Response Data.” Computer Methods and Programs in
Biomedicine, 58, 25–34.
Yee T (2010). “The VGAM Package for Categorical Data Analysis.” Journal of Statistical
Software, 32, 1–34. URL http://www.jstatsoft.org/v32/i10/.
Yu K, Yuan W (2004). “Regression Models for Unbalanced Longitudinal Ordinal Data: Computer Software and a Simulation Study.” Computer Methods and Programs in Biomedicine,75, 195–200.

Affiliation:
Anestis Touloumis
University of Brighton
School of Architecture, Technology and Engineering
University of Brighton
Brighton, BN2 4GJ, UK
E-mail: A.Touloumis@brighton.ac.uk

Journal of Statistical Software
published by the Foundation for Open Access Statistics
March 2015, Volume 64, Issue 8
doi:10.18637/jss.v064.i08

https://www.jstatsoft.org/
https://www.foastat.org/
Submitted: 2013-06-08
Accepted: 2014-07-24