5

Ok, I have searched and searched and just have no clue where to start. First, what I would like to do is produce a QQ-plot (or even a readable residual plot) to look at the fit of my model. I guess I just don't understand how the parameters that go into qnbinom() are obtained from the output of MASS::glm.nb(). I am attempting to use probplot() from package e1071, but am unsure of the inputs needed. It would be great if someone having experience with fitting negative binomials could lend a hand.

Secondly, I came across a residual plot here: http://www.stat.cmu.edu/~hseltman/Rclass/R8.R I can make it to work, but I don't know how to interpret it or if I am using it correctly. Has anyone else used this?

At the moment I am relying on the AIC and plots of fitted vs. actual values to assess the fit of my model and I would like something a little better!

Edit: Hopefully this will clarify what I am asking. With qnbinom(p, size, prob, mu, lower.tail = TRUE, log.p = FALSE), how do I (or is it even possible to) get the p, size, prob, mu from the output of a glm.nb fit model? From my research, I have found that size is the dispersion parameter, but other than that I'm not sure where to go. I know theta goes in there somehow, just not sure how to get it in the form needed.

Edit 2: Ok, once I have a distplot(), is there a guide to interpreting it? I am fairly positive I have a bad fit because I have a curved plot with a red line going through it (with many points on the tails far from the red line). The prob: ML = 0.011, is this rejecting that the distribution is from the negative binomial specified?

user1228982
  • 61
  • 1
  • 4
  • Regarding your edit: you only need size and prob OR size and mu. In the former case size is the parameter for target for number of successful trials, in the latter it is the dispersion parameter, namely theta! Both versions are related by prob = size/(size+mu). You can extract theta from the glmnbobject by object$theta and the mean you get by exp(coef(object)). See the example below. If you have a regression model, getting the mean is a bit trickier, you need to put in all the values for the x as well. Is the latter where the problems lie? See below for an example. – Momo Mar 28 '12 at 21:05
  • @Momo thanks for the detailed explanation. I do understand I have a model which requires transformation. This is a regression model with multiple predictors (some factors, others continuous predictors), so taking exp(coef()) just gives me a list. Would I just multiply those all together? Wish this distributuon was covered more thoroughly in one of my courses, but we (and our book) stopped with the Poisson and Binomial... So, I've been left scouring the internet for info. For the p in qnbinom(), what is supposed to be used? I don't understand what probabilities are needed. Thanks! – user1228982 Mar 29 '12 at 13:54
  • After a little research, do I get p by using ppoints(data)? (obtained from: http://stats.stackexchange.com/questions/4987/understanding-a-quantile-quantile-plot) – user1228982 Mar 29 '12 at 14:02
  • One more question: In the output of summary(model), it shows sometihng like the following: (Dispersion parameter for Negative Binomial(10.1883) family taken to be 1). I thought this meant the dispersion parameter was 1, while theta was 10.1883. You seemed to imply both are the same. Am I wrong in that interpratation? Also, I found this post useful, for anyone that might reference this in the future: http://intellinexus.wordpress.com/2010/11/29/creating-a-q-q-plot/ – user1228982 Mar 29 '12 at 14:45
  • You can find an explanation below what to do with getting the mean from the model with predictors, the second edit. But that will give you the same as fitting the intercept only model exp(coef(glm.nb(y~1))), which is the mean mu. p are the probabilities for which you want the quantiles which you must supply and does not come from glm.nb or nayone else. If you want to know the 95% Quantile, supply 0.95. – Momo Mar 29 '12 at 19:52
  • About the dispersion confusion: Given theta, the negative binomial distribution is an exponential family with a single parameter hence the dispersion parameter is set to 1 (as in the Poisson or Binomial). In the model with unknown theta, theta is also a dispersion parameter (called that way in R's nbinom) but you may call it shape parameter to avoid confusion. – Momo Mar 29 '12 at 20:03
  • let us [continue this discussion in chat](http://chat.stackexchange.com/rooms/2938/discussion-between-momo-and-user1228982) – Momo Mar 29 '12 at 20:05
  • Not enough reps to enter chat, but your comment cleared things up for me. I couldn't find anywhere that really clarified the parameter naming in R when using nbinom, which was causing some headaches. Thanks so much! – user1228982 Apr 06 '12 at 16:01
  • `glm.nb` provides the respective fitted means as: > unique(m1$fitted.values) [1] 17.95455 15.22500 – Mischief_Monkey Jun 28 '17 at 08:11

1 Answers1

6

You might find distplot() from the vcd package useful either for the original data (edit: you can't use it on residuals). This plots Friendly's "negativebinomialness plots" and provides how well the negative binomial model fits
distplot(response, type = "nbinomial", ...)

To obtain the parameters: glm.nb uses the "Gamma mixture of Poisson" representation. It is actually a log-linear model that is fitted, so you should get the mean as $\exp(X\beta)$.

For example, let's say your data come from a negbin with mean 5 and theta of 1 (in the alternative representation as described above). Then you can get the mean estimate simply by 

set.seed(10)  
df <- data.frame(y=rnbinom(100,size=1,mu=5))  
m0 <- glm.nb(y~1,data=df)  
m0  
exp(coef(m0))  
m0$theta

which are in this case 5.1 for the mean (pretty close) and 1.6 for the dispersion parameter (pretty far off).

If you fit a model for the conditional mode, you interpret it accordingly as in every other log linear model, see this discussion on stack exchange.

EDIT: If you want to know how to get the mean in a negbin regression model you need to sum up the linear predictor $X\beta$.

For example: I take the quine data and fit 

m1 <- glm.nb(Days~Sex,data=quine)

now males are 1 females are 0. To get the mean for males you write 

> exp(coef(m1)[1]+coef(m1)[2]*1)  
[1] 17.95455

and for females

> exp(coef(m1)[1]+coef(m1)[2]*0)     
[1] 15.225

Now to get the mean you must weight this with the occurence of all females and males which is 

> table(quine$Sex)  
 F  M   
80 66

and hence the mean is

> (80/(66+80))*15.225+(66/(80+66))*17.95455  
[1] 16.45685

This is confirmed by

> nb0 <- glm.nb(Days ~ 1, data = quine)    
> exp(coef(nb0))  
(Intercept)  
[1] 16.4589

(apart from rounding errors).

Community
  • 1
Momo
  • 8,839
  • 3
  • 46
  • 59
  • Thanks Momo, I understand that the negative binomial is a mixture, it's just using the numbers I am given in the glm.nb output that is giving me trouble. I saw the post you mentioned, but didn't see what I was looking for. distplot looks promising, however when I tried to use it I get an error which I can't make sense of: Error in nfreq[count + 1] – user1228982 Mar 28 '12 at 19:12
  • What exactly is the problem in the output? I don't quite know what it is you need to know. You get a linear predictor which you need to to take exp($X\beta$) because your model is $log(\mu)=X\beta$ with mu being the expected count. So in fact you fit this model: $mu=exp(\beta_0)*exp(\beta_1x_1)...$. This is how to interpret them, as a loglinear model. Regarding distplot, I don't know the error, but perhaps you are right it might only take the dependent variable as an input (which must be 0 or positive integer), hence residuals may be not allowed. Can't you just use the dep. variable though? – Momo Mar 28 '12 at 20:22
  • I think my main problem is how the parameters are named in different functions vs the output from glm.nb(). – user1228982 Mar 29 '12 at 15:31
  • Like I said: mu is exp(coef(glm.nb(y~1))) with y being the dependent variable. size is theta. prob is (size/size+mu). – Momo Mar 29 '12 at 19:46
  • 1
    Ok, so I could really just use the mean() function on my original data for mu... that's what I was thinking but for some reason it seemed too simple. I think I skipped over this line: nb0 – user1228982 Mar 29 '12 at 20:25