# R/Stata package for zero-truncated negative binomial GEE?

this is my first post. I'm truly grateful for this community.

I am trying to analyze longitudinal count data that is zero-truncated (probability that response variable = 0 is 0), and the mean != variance, so a negative binomial distribution was chosen over a poisson.

Functions/commands I've ruled out:

R

• gee() function in R does not account for zero-truncation nor the negative binomial distribution (not even with the MASS package loaded)
• glm.nb() in R doesn't allow for different correlation structures
• vglm() from the VGAM package can make use of the posnegbinomial family, but it has the same problem as Stata's ztnb command (see below) in that I can't refit the models using a non-independent correlation structure.

Stata

• If the data wasn't longitudinal, I could just use the Stata packages ztnb to run my analysis, BUT that command assumes that my observations are independent.

I've also ruled out GLMM for various methodological/philosophical reasons.

For now, I've settled on Stata's xtgee command (yes, I know that xtnbreg also does the same thing) that takes into account both the nonindependent correlation structures and the neg binomial family, but not the zero-truncation. The added benefit of using xtgee is that I can also calculate qic values (using the qic command) to determine the best fitting correlation structures for my response variables.

If there is a package/command in R or Stata that can take 1) nbinomial family, 2) GEE and 3) zero-truncation into account, I'd be dying to know.

I'd greatly appreciate any ideas you may have. Thank you.

-Casey

For R two options spring to mind, both of which I am only vaguely familiar with at best.

The first is the pscl package, which can fit zero truncated inflated and hurdle models in a very nice, flexible manner. The pscl package suggests the use of the sandwich package which provides "Model-robust standard error estimators for cross-sectional, time series and longitudinal data". So you could fit your count model and then use the sandwich package to estimate an appropriate covariance matrix for the residuals taking into account the longitudinal nature of the data.

The second option might be to look the geepack package which looks like it can do what you want but only for a negative binomial model with known theta, as it will fit any type of GLM that R's glm() function can (so use the family function from MASS).

A third option has raised it's head: gamlss and it's add-on package gamlss.tr. The latter includes a function gen.trun() that can turn any of the distributions supported by gamlss() into a truncated distribution in a flexible way - you can specify left truncated at 0 negative binomial distribution for example. gamlss() itself includes support for random effects which should take care of the longitudinal nature of the data. It isn't immediately clear however if you have to use at least one smooth function of a covariate in the model or can just model everything as linear functions like in a GLM.

• The pscl package, I believe, only fits zero-inflated and hurdle models. Hurdle models incorporate both a left-truncated count component and a right-censored hurdle component. I don't how or even if I am able to run a hurdle model without the hurdle component, but I will look into the sandwick package. As for the geepack package, it seems to have the same problem as gee package; when I specify a "negative.binomial" family (from MASS), without specifying a theta, it will ask for a theta. However, when I specify a theta value, it will spit out an error saying it's an unrecognized family. – Iris Tsui Feb 20 '11 at 15:52
• @Casey - sorry I misread your requirements re zero truncation. Shame that geepack doesn't work with that family function. If I think of anything else, I'll update here. – Gavin Simpson Feb 20 '11 at 16:56
• @Casey I've added a note about the gamlss package which might fit the bill in R also. – Gavin Simpson Feb 21 '11 at 10:27
• Accepting your answer because of the multiple suggestions for resources and functions that have improved my understanding. It seems like 'gamlss' would be a possible way to solve my issue, but because I'm actually a non-statistician, I don't currently have the background in math nor the time to open that can of worms right now (but maybe eventually I will). As mentioned in another comment, for my data at least, it seems that ignoring zero-truncation won't change my estimates and std errors much. For my intended audience, I believe a nbinomial GEE will do just fine. Thanks! – Iris Tsui Feb 25 '11 at 0:32

Hmm, good first question! I don't know of a package that meets your precise requirements. I think Stata's xtgee is a good choice if you also specify the vce(robust) option to give Huber-White standard errors, or vce(bootstrap) if that's practical. Either of these options will ensure the standard errors are consistently estimated despite the model misspecification that you'll have by ignoring the zero truncation.

That leaves the question of what effect ignoring the zero truncation will have on the point estimate(s) of interest to you. It's worth a quick search to see if there is relevant literature on this in general, i.e. not necessarily in a GEE context -- i would have thought you can pretty safely assume any such results will be relevant in the GEE case too. If you can't find anything, you could always simulate data with zero-truncation and known effect estimates and assess the bias by simulation.

• I made sure to estimate robust standard errors. Also, in the book "Mixed effects models and extensions in ecology with R" by Zuur, et al, 2009, on page 261, they mention, "if the mean of the response variable is relatively large, ignoring the truncation problem, then applying a Poisson or negative binomial (NB) generalised linear model (GLM), is unlikely to cause a problem." Fortunately, the means of my response variables are large, so I feel a little more comfortable deprioritizing zero-truncation compared to the GEE and negbinomial aspects of my regressions. – Iris Tsui Feb 20 '11 at 7:28
• Sounds like you already know more about this topic than I do! Or anyone else on this site, judging by the lack of other responses. – onestop Feb 20 '11 at 8:29
• It is a little unbelievable; who knew that overdispersed longitudinal count data would be so difficult to analyze (without doing a GLMM, which I haven't even looked into doing yet)? If only my data were zero-inflated, that would be another story. – Iris Tsui Feb 20 '11 at 8:46

I had the same issue in my dissertation. In Stata, I just built myself a custom .ado program with two calls to xtgee.

For this, I found the "Modeling Health Care Costs and Counts" slides/programs by Partha Deb, Willard Manning, and Edward Norton to be useful. They don't talk about longitudinal data, but it's a useful starting point.

I was looking for answers on glmmADMB interpretation and I saw your post. I know it was a long time ago but I might have the answer.

Look into the package glmmADMB when using hurdle models. You have to split into two the analyses of your data: one of them treats just the no zero data. You may add mixed effects and chose the distribution. The condition is that the data has to be zero-inflated and I do not know if this fitted your requirements! Anyways, I hope you found out long ago!