Any updates to model negative binomial distribution data with GEE in R?

Question

Goal: run GEE with negative binomial distribution to analyze dataset in R

This question was asked here, but the answers are ~5 years old and I wonder if there are new developments. I summarize my data and then outline how I've tried to find a solution.

My model (successfully runs using a Poisson rather than negative binomial distribution with package geeglm):

m1 <- geeglm(Diract ~ Dir*Rec + Year + offset(LnScan), 
             family = poisson("log"), data = Direct, id = ID, corstr = "exchangeable")

Data are here

Diract: interval outcome variable, counts of behaviors from directors (“Dir” in group B or C) to recipients (“Rec" of groups A, B, or C). Range from 0-4, overdispersed.
Dir*Rec: main predictor of interest, interaction of two factors (director and recipient group). Animals belong to one of 3 groups (life stages): A, B, C.
Year: year of data collection (2010, 2011, 2012) Goal family: negative binomial with log link function

Some individuals (~5) transitioned from one group to another between periods of data collection during the three years of study. I have multiple measures for each director (both within and between years). E.g., an individual who was alive for the entire study has counts of the behavior directed toward groups A, B, and C for each of the three years of study (total = 9 counts). Some individuals were observed all three years, others were only observed one or two years.

Initially I used GEE in SPSS, software I longer have access to. In SPSS, QIC values indicated that negative binomial was a dramatically better fit for my data vs. Poisson.

Solution1 suggested via Stackoverflow link above: use library("sos") and findFn("{generalized estimating equation}") to find packages in R that run GEEs. I checked all 14 packages suggested. geepack lacks negative binomial family; Zelig, JGEE, pseudo, aftgee, etm, wgeesel, MethylCapSig, and miLineage leverage geepack. MESS, PGEE, spind, and threeboost lack negative binomial, and repolr is solely for ordinal outcomes.

Solution2 suggested by Ben Bolker on Nabble: "I would try glmmPQL in the MASS package. I don't think you can quite get negative binomial regression this way, but you can definitely get a quasipoisson model. I think exchangeable correlation corresponds to correlation=corCompSymm() in your glmmPQL command."

My attempt with error message:

m2 <- glmmPQL(Diract ~ Dir*Rec + offset(LnScan) + Year, 
              random = ~ 1 | ID, 
              family = negative.binomial(1), 
              data = Direct, 
              correlation = corCompSymm(form=~1|ID))

Error in glmmPQL(Diract ~ DirPar * RecPar + offset(LnScan) + Year, random = ~1 | : could not find function "corCompSymm"

If anyone has tips or solutions, I'm open to anything. Thank you very much!

hghg hghg hghg hghg · Accepted Answer · 2019-07-26T06:59:39

This question has been asked a number of times over the years and I came across it after encountering the same problem and searching for an answer for a couple of hours. Two solutions:

The geeM package allows use of the negative binomial family
via the R reticulate package import the python statsmodels module and call the gee fitting function within R which allows the negative binomial family

Any updates to model negative binomial distribution data with GEE in R?

1 Answers