[R] negative binomial family glm R and STATA

Mikis Stasinopoulos d.stasinopoulos at londonmet.ac.uk
Mon Jan 8 00:14:05 CET 2007

Dear  Patrick

Try the package gamlss which allow to fit the negative binomial 
distrbution. For example  with your data I am getting
GAMLSS-RS iteration 1: Global Deviance = 817.9027
GAMLSS-RS iteration 2: Global Deviance = 817.9025

Family:  c("NBI", "Negative Binomial type I")
Fitting method: RS()

Call:  gamlss(formula = nbcas ~ ., family = NBI, data = zonesdb4)

Mu Coefficients:
 (Intercept)           pop          Area    V_100kHab1        gares1 
   3.204e+00     1.114e-05     1.354e-05     9.144e-01     7.946e-01 
      ports1  axe_routier1  axe_routier2         lacs1 
  -1.730e+00     1.989e-01            NA     3.042e+00 
Sigma Coefficients:

 Degrees of Freedom for the fit: 9 Residual Deg. of Freedom   83
Global Deviance:     817.902
            AIC:     835.902
            SBC:     858.599

Note that the AIC: 835.902 is similar to your fitted model using  glm.nb 
which is  AIC: 836.2.

The coefficients are not identical but this is not suprissing when you 
are using x-variables with extreme values as pop and Area.
The profile function for sigma can be found using

prof.dev(ga1,"sigma", min=7, max=16, step=0.1, type="l")

Your discrepancy with STATA come from the fact that in STATA you are 
fitting the model with sigma fixed to 1.
You can see that by fitting the same model in  GAMLSS.

 >  ga2<-gamlss(nbcas~.,data=zonesdb4,family=NBI, sigma.fix=T, 
GAMLSS-RS iteration 1: Global Deviance = 1194.299
GAMLSS-RS iteration 2: Global Deviance = 1194.298

This  is similar  to  the  log likelihod you are  getting in STATA. i.e. 
-2*-597.14778= 1194.296.

You can also use the stepGAIC() function to simplify your model. For  


will result to a model with only pop and lacs in the mdel.

Note also the the Negative binomial type II fits better to you data.

 >  ga3<-gamlss(nbcas~.,data=zonesdb4,family=NBII)
GAMLSS-RS iteration 1: Global Deviance = 804.5682
GAMLSS-RS iteration 10: Global Deviance = 804.4995



More information about the R-help mailing list