[R] negative binomial family glm R and STATA

From: Patrick Giraudoux <patrick.giraudoux_at_univ-fcomte.fr>
Date: Sat 06 Jan 2007 - 11:00:26 GMT


Dear Lister,

I am facing a strange problem fitting a GLM of the negative binomial family. Actually, I tried to estimate theta (the scale parameter) through glm.nb from MASS and could get convergence only relaxing the convergence tolerance to 1e-3. With warning messages:

 glm1<-glm.nb(nbcas~.,data=zonesdb4,control=glm.control(epsilon = 1e-3)) There were 25 warnings (use warnings() to see them)
> warnings()

Warning messages:
1: iteration limit reached in: theta.ml(Y, mu, n, w, limit = control$maxit, trace = control$trace > ... 2: NaNs produced in: sqrt(1/i)

etc....

The estimate of theta was: 0.0939. When trying to compute confidence interval then, I got this message:

> confint(glm1a)

Waiting for profiling to be done...
Error in profile.glm(object, which = parm, alpha = (1 - level)/4, trace = trace) :

        profiling has found a better solution, so original fit had not converged

Moreover, trying some other solutions "by hand" (without warnings produced, here) with glm(.... family=negative.binomial(1)....), I found that theta = 0.7 lead to a much lower AIC (AIC= 1073) than theta = 1 (AIC=1211). Facing such unstable results my first reaction has been to forget about fitting a negative binomial model on this data set. The reader will find the dataset in a dumped format below for trials.

A friend of mine tried the same with STATA and got the following result without any warning from STATA :

. glm nbcas pop area v_100khab gares ports axe_routier lacs, family(nbinomial) link(log) eform  

Iteration 0:   log likelihood = -616.84031 
Iteration 1:   log likelihood = -599.77767 
Iteration 2:   log likelihood = -597.22486 
Iteration 3:   log likelihood = -597.14784 
Iteration 4:   log likelihood = -597.14778 
Iteration 5:   log likelihood = -597.14778 
 
Generalized linear models                          No. of obs      
=        92
Optimization     : ML                              Residual df     
=        84
                                                   Scale parameter 
=         1
Deviance         =  597.0007978                    (1/df) Deviance =  
7.107152
Pearson          =  335.6135273                    (1/df) Pearson  =  
3.995399  
Variance function: V(u) = u+(1)u^2                 [Neg. Binomial]
Link function    : g(u) = ln(u)                    [Log]
 
                                                   AIC             =  
13.15539
Log likelihood   = -597.1477759                    BIC             =  
217.1706  
             |                 OIM
       nbcas |        IRR   Std. Err.      z    P>|z|     [95% Conf. 
Interval]
-------------+---------------------------------------------------------------- 

         pop |   1.000011   1.82e-06     6.02   0.000     1.000007    
1.000014
        area |   1.000014   .0000244     0.57   0.569     .9999661    
1.000062
   v_100khab |   2.485394   .7924087     2.86   0.004     1.330485    
4.642806
       gares |   2.185483   .7648255     2.23   0.025     1.100686    
4.339418
       ports |   .1805793    .100423    -3.08   0.002     .0607158    
.5370744
 axe_routier |    .828243   .2258397    -0.69   0.489     .4853532    
1.413376
        lacs |   20.50183   12.17126     5.09   0.000     6.404161    
65.63311

Has somebody an idea about (1) why the AIC values given are so different between softwares (R = 1211, STATA= 13.15) for the same model and (2) what can explain so different behaviour between R and STATA ?

Here below the data.frame:

zonesdb4 <-
structure(list(nbcas = as.integer(c(318, 0, 42, 3011, 6, 911,

45, 273, 0, 0, 89, 122, 407, 83, 0, 1844, 58, 0, 0, 0, 0, 8926,
0, 0, 0, 0, 108, 0, 13, 1884, 0, 0, 0, 0, 963, 0, 199, 735, 0,
2182, 971, 0, 65, 0, 7927, 30, 0, 186, 0, 1363, 808, 0, 0, 0,
0, 135, 0, 1338, 0, 0, 488, 0, 344, 0, 0, 454, 4808, 0, 692,
0, 0, 24, 1301, 0, 0, 474, 228, 0, 0, 98, 44, 0, 0, 0, 1562,
375, 327, 0, 0, 0, 0, 0)), pop = as.integer(c(247215, 55709, 63625, 253153, 51789, 142806, 129839, 95799, 129996, 66668, 76043, 267232, 153200, 136333, 264888, 198244, 233600, 89152, 128085, 71803, 81911, 122523, 149806, 122470, 50979, 160773, 80700, 56146,
226965, 245322, 165768, 215129, 46843, 108471, 108690, 188724,
161794, 226965, 95850, 156326, 145291, 51789, 218808, 53189,
245854, 152047, 146509, 243765, 65012, 226830, 66742, 144762,
93858, 73793, 123107, 189793, 91013, 135212, 67487, 105050, 194903, 206606, 62169, 96832, 145570, 167062, 1598576, 146509, 103928,
118334, 91509, 295644, 139650, 106980, 66529, 126126, 257341,
56973, 33793, 164072, 145225, 155638, 131100, 100880, 245482,
166213, 127365, 113713, 57540, 78571, 62499, 81916)), Area = c(10027.1,
9525.3, 638.646, 861.486, 4966.32, 11973, 1823.89, 1327.45, 789.595,
4892.38, 638.908, 15959.8, 2036.56, 7397.62, 4626.03, 10237.5,
9823.24, 4253.59, 2448.78, 12280.2, 910.972, 16365, 2041.92,
4343.46, 1081.42, 1601.11, 4664.47, 335.865, 2818.68, 7348.1,
1148.41, 265.158, 14883.6, 3698.58, 12444.6, 1711.45, 15462,
10419.5, 13119.2, 1276.14, 872.91, 19291.4, 6719.82, 8505.53,
13219.8, 13069, 5212.03, 3924.42, 791.219, 881.281, 10038.5,
9101.94, 7925.71, 943.062, 9888, 20585.3, 4600.35, 3258.27, 11813.4,
130.184, 10644.3, 1925.89, 1892.88, 3833.6, 350.3, 7154.79, 2800.63,
559.986, 3152, 7095.39, 6058.3, 113.225, 5067.66, 1293.05, 15109.8,
4111.54, 94.5213, 4012.91, 1468.02, 10651.3, 8541.69, 1806.28,
20166.3, 1110.75, 2026.98, 21114.4, 2041.51, 17740.9, 16627.5,
15266.1, 525.467, 371.132), V_100kHab = structure(as.integer(c(1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2,
1, 1, 1, 1, 1, 1, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 1, 2,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 1, 1, 1, 1, 2, 1, 1, 1, 1,
1, 1, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 2, 1, 1, 1, 1, 2)), .Label = c("0", "1"), class = "factor"),

    gares = structure(as.integer(c(2, 2, 1, 1, 1, 1, 1, 1, 1,

    2, 2, 2, 1, 1, 1, 2, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 2, 1,
    1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
    1, 2, 1, 2, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1,
    1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 2, 1, 2, 2, 1, 1,
    1, 1, 1, 1, 1, 1, 1)), .Label = c("0", "1"), class = "factor"),
    ports = structure(as.integer(c(2, 1, 1, 1, 1, 1, 1, 1, 1,
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1,
    1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2,
    2, 1, 1, 1, 1, 2, 1, 2, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1,
    1, 1, 2, 1, 1, 1, 1, 1, 1, 2, 1, 2, 1, 1, 2, 1, 2, 1, 1,
    2, 1, 1, 1, 1, 1, 1)), .Label = c("0", "1"), class = "factor"),
    axe_routier = structure(as.integer(c(2, 3, 3, 3, 2, 2, 2,
    3, 2, 3, 2, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 3, 2, 2, 3, 2,
    3, 2, 3, 3, 3, 2, 2, 3, 2, 3, 2, 3, 2, 3, 3, 2, 3, 2, 3,
    2, 2, 3, 3, 2, 3, 2, 2, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2,
    3, 3, 2, 2, 3, 3, 3, 3, 2, 3, 3, 2, 3, 3, 2, 3, 2, 3, 3,
    3, 2, 3, 3, 3, 2, 2, 3, 3)), .Label = c("0", "1", "2"), class = 
"factor"),
    lacs = structure(as.integer(c(1, 1, 1, 1, 1, 1, 1, 1, 1,
    1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1,
    1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 2,
    2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1,
    1, 1, 2, 1, 1, 1, 1, 1, 1, 2, 1, 2, 1, 1, 1, 1, 1, 1, 1,
    2, 1, 1, 1, 1, 1, 1)), .Label = c("0", "1"), class = "factor")), .Names = c("nbcas",
"pop", "Area", "V_100kHab", "gares", "ports", "axe_routier",
"lacs"), row.names = c("Ankoro", "Dilolo", "Tshitenge", "Tshitshimbi",
"Tshofa", "Tshudi Loto", "Vanga Kete", "Wikong", "Djalo Djeka",
"Fungurume", "Gandajika", "Kabalo", "Kabeya Kamuanga", "Kabinda",
"Kabondo Dianda", "Kabongo", "Kafakumba", "Bena Dibele", "Kafubu",
"Kalamba", "Kalambayi Kabanga", "Kalemie", "Kalenda", "Kalonda Est",
"Kamalondo", "Kamana", "Kambove", "Kamiji", "Bibanga", "Kamina",
"Kampemba", "Kanda Kanda", "Kansimba", "Kanzenze", "Kapanga",
"Kapolowe", "Kasaji", "Kasenga", "Katako Kombe", "Katuba", "Kenya",
"Kiambi", "Kikula", "Kilela Balanda", "Kilwa", "Kinda", "Bukenya",
"Kinkondja", "Kipushi", "Kisanga", "Kitenge", "Kole", "Kongolo",
"Likasi", "Lodja", "Lomela", "Lualaba", "Butumba", "Lubao", "Lubilanji",
"Lubudi", "Lubumbashi", "Ludimbi Lukula", "Lukafu", "Lukelenge",
"Lusambo", "Cilundu", "Makota", "Malemba Nkulu", "Manika", "Manono",
"Miabi", "Minga", "Muene Ditu", "Mufunga Sampwe", "Mukanga",
"Mukumbi", "Mulongo", "Mulumba", "Mutshatsha", "Nyemba", "Dilala",
"Nyunzu", "Panda", "Pania Mutombo", "Pweto", "Rwashi", "Sakania",
"Sandoa", "Songa", "Tshamilemba", "Tshilenge"), class = "data.frame", 
na.action = structure(as.integer(c(8,
34, 40, 41, 45, 71, 73, 79, 80, 83, 84, 85, 86, 93)), .Names = c("Wembo Nyama",
"Kaniama", "Kasansa", "Bukama", "Kayamba", "Luputa", "Lwamba",
"Mbuji mayi", "Mbulula", "Mitwaba", "Moba", "Dikungu Tshumbe",
"Mpokolo", "Mumbunda"), class = "omit"))


	[[alternative HTML version deleted]]

______________________________________________
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 Sat Jan 06 22:10:36 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 Sun 07 Jan 2007 - 19:31:41 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.