From: Mikis Stasinopoulos <d.stasinopoulos_at_londonmet.ac.uk>

Date: Sun 07 Jan 2007 - 23:14:05 GMT

(Intercept)

R-help@stat.math.ethz.ch mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Mon Jan 08 17:34:55 2007

Date: Sun 07 Jan 2007 - 23:14:05 GMT

Dear Patrick

Try the package gamlss which allow to fit the negative binomial distrbution. For example with your data I am getting

#---------------------------------------------------ga1<-gamlss(nbcas~.,data=zonesdb4,family=NBI) GAMLSS-RS iteration 1: Global Deviance = 817.9027 GAMLSS-RS iteration 2: Global Deviance = 817.9025 ga1

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+00Sigma Coefficients:

(Intercept)

2.313

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,

sigma.start=1)

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 example

ga2<-stepGAIC(ga1)

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

Thanks

Mikis

R-help@stat.math.ethz.ch mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Mon Jan 08 17:34:55 2007

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.1.8, at Mon 08 Jan 2007 - 08:30:26 GMT.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help.
Please read the posting
guide before posting to the list.
*