# Multinomial logistic regression with geepack in R

I am working on fitting a GEE model to a multinomial logistic outcome using the R package geepack. My understanding is that the package uses glm to fit the formula. However, from what I can see, glm doesn't have a family for a multinomial outcome, just binomial. How do you go about fitting a multinomial logistic equation using gee in R? Thanks.

• This question seems to be only about how to do this in R. Thus it is off-topic for CV (see our help center). If you have a question about the statistical aspects involved, please edit to clarify. Otherwise, this Q may be closed. – gung - Reinstate Monica Aug 1 '14 at 18:03
• This question appears to be off-topic because it is about how to use R. – gung - Reinstate Monica Aug 1 '14 at 18:04
• @JessicaW. Dont's use glmnet unless you have no other option. This package is not particularly user friendly and intuitive. Your most user friendly option is package nnet. Here is a very well exaplained example of using the multinom function ats.ucla.edu/stat/r/dae/mlogit.htm Also, here is my question about MLR using R in case it is helpful stackoverflow.com/questions/22293517/… – Koba Aug 1 '14 at 18:20

Try using the R package glmnet. It has very good documentation, it is the state of the art and certainly does multinomial regression. I have noticed some of these packages expect you to figure out cross validation yourself. In that sense, glmnet is nice. The cv function does it for you.

Besides I don't know what algorithm GEE uses (I haven't looked into it, they probably do the newer methods as is). Some older packages resort to Interior Point Methods which was state of the art at the time but has since become outdated.

• Would I use this to set up the model alone, and then use geepack to run the generalized estimating equation? Or would glmnet also run the generalized estimating equation? – Jess A. Aug 1 '14 at 18:51

Another option is the multgee package. The included example of a nominal (there's also an ordinal example) multinomial logistic GEE-solved marginal regression model predicts housing status (y=1 for "street living", 2 for community living, and 3 for independent housing) from time in months and Section 8 rent certificate status (sec, binary). Here's the code: require(multgee);data(housing) summary(nomLORgee(y~factor(time)*sec,data=housing,id=id,repeated=time)) And selected output:

Link : Baseline Category Logit

Local Odds Ratios:
Structure:         time.exch
Model:             3way
Homogenous scores: TRUE

Summary of residuals:
Min.   1st Qu.    Median      Mean   3rd Qu.      Max.
-0.585900 -0.278000 -0.128700 -0.002166  0.421200  0.908800

Number of Iterations: 3

Coefficients:
Estimate   san.se   san.z Pr(>|san.z|)
beta01                1.66073  0.25026  6.6359      < 2e-16 ***
factor(time)6:1      -1.87010  0.31876 -5.8668      < 2e-16 ***
factor(time)12:1     -2.92506  0.36829 -7.9424      < 2e-16 ***
factor(time)24:1     -2.81359  0.34258 -8.2130      < 2e-16 ***
sec:1                -0.53680  0.33704 -1.5927      0.11122
factor(time)6:sec:1  -1.18218  0.46036 -2.5680      0.01023 *
factor(time)12:sec:1  0.07916  0.48306  0.1639      0.86984
factor(time)24:sec:1  0.03273  0.46558  0.0703      0.94396
beta02                1.16643  0.26273  4.4397        1e-05 ***
factor(time)6:2      -0.25454  0.30080 -0.8462      0.39743
factor(time)12:2     -0.57052  0.31176 -1.8300      0.06725 .
factor(time)24:2     -1.04101  0.30716 -3.3892      0.00070 ***
sec:2                -0.10704  0.34759 -0.3080      0.75811
factor(time)6:sec:2  -1.62342  0.41349 -3.9261        9e-05 ***
factor(time)12:sec:2 -2.04850  0.44543 -4.5990      < 2e-16 ***
factor(time)24:sec:2 -1.04965  0.41831 -2.5093      0.01210 *

Local Odds Ratios Estimates:
[,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]
[1,] 0.000 0.000 1.359 1.882 1.359 1.882 1.359 1.882
[2,] 0.000 0.000 1.882 3.687 1.882 3.687 1.882 3.687
[3,] 1.359 1.882 0.000 0.000 1.359 1.882 1.359 1.882
[4,] 1.882 3.687 0.000 0.000 1.882 3.687 1.882 3.687
[5,] 1.359 1.882 1.359 1.882 0.000 0.000 1.359 1.882
[6,] 1.882 3.687 1.882 3.687 0.000 0.000 1.882 3.687
[7,] 1.359 1.882 1.359 1.882 1.359 1.882 0.000 0.000
[8,] 1.882 3.687 1.882 3.687 1.882 3.687 0.000 0.000

pvalue of Null model: <0.0001


See the following reference for background and interpretation of this model. Also, here's the structure of the housing dataset (using str(lapply(housing,factor))) in case it helps you organize yours:

List of 4
$id : Factor w/ 362 levels "1","2","3","4",..: 1 1 1 1 2 2 2 2 3 3 ...$ y   : Factor w/ 3 levels "0","1","2": 2 3 3 3 2 3 3 2 1 3 ...
$time: Factor w/ 4 levels "0","6","12","24": 1 2 3 4 1 2 3 4 1 2 ...$ sec : Factor w/ 2 levels "0","1": 2 2 2 2 2 2 2 2 2 2 ...


Reference
Touloumis, A. (2011). Generalized estimating equations for multinomial responses. PhD dissertation, University of Florida. Retrieved from http://ufdcimages.uflib.ufl.edu/UF/E0/04/32/26/00001/touloumis_a.pdf
.

• Hello, I read the example on p.182-185 of the thesis but didn't quite get it... Could anyone help a little more about the interpretation? Thank you! – leoce Sep 12 '18 at 4:07