Running 'test_irls.R' [152s]
  Running 'test_ml_g.R' [161s]
  Running 'test_ml_glm.R' [1s]
  Running 'test_ml_glm2.R' [159s]
  Running 'test_ml_glm3.R' [1s]
  Running 'test_nbinomial.R' [174s]
Running the tests in 'tests/test_irls.R' failed.
Complete output:
  > library(msme)
  Loading required package: MASS
  Loading required package: lattice
  > 
  > # library(msme, lib.loc="lib")
  > 
  > library(MASS)
  > 
  > data(medpar)
  > 
  > denom <- rep(1:5, each=299, times=1)*100   # m : binomial denominator w medpar
  > oset <- rep(1:5, each=299, times=1)*100    # offset Poisson, NB, offset w medpar
  > loset <- log(oset)                         #    log of oset
  > 
  > ## POISSON -------------------------------------------------
  > 
  > irls.poi <- irls(los ~ hmo + white,
  +                  family = "poisson",
  +                  link = "log",
  +                  data = medpar)
  > 
  > glm.poi <- glm(los ~ hmo + white,
  +                family = "poisson",
  +                data = medpar)
  > 
  > irls.poi
  
  Call:  irls(formula = los ~ hmo + white, data = medpar, family = "poisson", 
      link = "log")
  
  Coefficients:
  X(Intercept)          Xhmo        Xwhite  
        2.4823       -0.1416       -0.1909  
  
  Degrees of Freedom: 1494 Total (i.e. Null);  1492 Residual
  Null Deviance:	    NA 
  Residual Deviance: 8813 	AIC: 14530
  > glm.poi
  
  Call:  glm(formula = los ~ hmo + white, family = "poisson", data = medpar)
  
  Coefficients:
  (Intercept)          hmo        white  
       2.4823      -0.1416      -0.1909  
  
  Degrees of Freedom: 1494 Total (i.e. Null);  1492 Residual
  Null Deviance:	    8901 
  Residual Deviance: 8813 	AIC: 14530
  > 
  > summary(irls.poi)
  
  Call:
  irls(formula = los ~ hmo + white, data = medpar, family = "poisson", 
      link = "log")
  
  Deviance Residuals:
     Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  -4.1197 -1.7487 -0.6213 -0.3073  0.9434 18.9477 
  
  Coefficients (all in linear predictor):
              Estimate     SE     Z        p    LCL     UCL
  (Intercept)    2.482 0.0259 95.89        0  2.432  2.5330
  hmo           -0.142 0.0237 -5.97 2.43e-09 -0.188 -0.0951
  white         -0.191 0.0273 -7.00  2.6e-12 -0.244 -0.1374
  
  Null deviance: NA  on  1494 d.f.
  Residual deviance: 8812.942  on  1492 d.f.
  Null Pearson:  on  1494 d.f.
  Residual Pearson: 11547.37  on  1492 d.f.
  Dispersion:
  AIC:  14534.09
  
  Number of optimizer iterations:  4 
  
  > summary(glm.poi)
  
  Call:
  glm(formula = los ~ hmo + white, family = "poisson", data = medpar)
  
  Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
  (Intercept)  2.48225    0.02589  95.890  < 2e-16 ***
  hmo         -0.14158    0.02373  -5.966 2.43e-09 ***
  white       -0.19089    0.02728  -6.998 2.60e-12 ***
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  (Dispersion parameter for poisson family taken to be 1)
  
      Null deviance: 8901.1  on 1494  degrees of freedom
  Residual deviance: 8812.9  on 1492  degrees of freedom
  AIC: 14534
  
  Number of Fisher Scoring iterations: 5
  
  > 
  > # Null deviance
  > 
  > irls(los ~ 1,
  +      family = "poisson",
  +      link = "log",
  +      data = medpar)$deviance
  [1] 8901.134
  > 
  > # cf
  > summary(glm.poi)$null.deviance
  [1] 8901.134
  > 
  > ## RATE POISSON
  > 
  > irls.rpoi <- irls(los ~ hmo + white,
  +                   family = "poisson",
  +                   link = "log",
  +                   offset = loset,
  +                   data = medpar,
  +                   verbose = 1)
  1 -2.99395 0.03960149 -0.3138124 3252.482 
  2 -3.140051 0.004622379 -0.2930442 -314.8172 
  3 -3.151896 -0.001205592 -0.2897326 -1.476491 
  4 -3.151968 -0.00127432 -0.2896969 -4.615321e-05 
  5 -3.151968 -0.001274325 -0.2896969 5.456968e-12 
  > glm.rpoi <- glm(los ~ hmo + white + offset(loset),
  +                 family = poisson,
  +                 data = medpar)
  > 
  > irls.rpoi
  
  Call:  irls(formula = los ~ hmo + white, data = medpar, family = "poisson", 
      link = "log", offset = loset, verbose = 1)
  
  Coefficients:
  X(Intercept)          Xhmo        Xwhite  
     -3.151968     -0.001274     -0.289697  
  
  Degrees of Freedom: 1494 Total (i.e. Null);  1492 Residual
  Null Deviance:	    NA 
  Residual Deviance: 11840 	AIC: 17560
  > glm.rpoi
  
  Call:  glm(formula = los ~ hmo + white + offset(loset), family = poisson, 
      data = medpar)
  
  Coefficients:
  (Intercept)          hmo        white  
    -3.151968    -0.001274    -0.289697  
  
  Degrees of Freedom: 1494 Total (i.e. Null);  1492 Residual
  Null Deviance:	    11940 
  Residual Deviance: 11840 	AIC: 17560
  > 
  > summary(irls.rpoi)
  
  Call:
  irls(formula = los ~ hmo + white, data = medpar, family = "poisson", 
      link = "log", offset = loset, verbose = 1)
  
  Deviance Residuals:
      Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
  -5.41202 -2.04944 -0.51239 -0.05207  1.60773 16.10967 
  
  Coefficients (all in linear predictor):
              Estimate     SE         Z        p     LCL     UCL
  (Intercept) -3.15197 0.0259 -121.7894        0 -3.2027 -3.1012
  hmo         -0.00127 0.0237   -0.0537    0.957 -0.0478  0.0452
  white       -0.28970 0.0273  -10.6179 2.46e-26 -0.3432 -0.2362
  
  Null deviance: NA  on  1494 d.f.
  Residual deviance: 11837.32  on  1492 d.f.
  Null Pearson:  on  1494 d.f.
  Residual Pearson: 16654.66  on  1492 d.f.
  Dispersion:
  AIC:  17558.47
  
  Number of optimizer iterations:  5 
  
  > summary(glm.rpoi)
  
  Call:
  glm(formula = los ~ hmo + white + offset(loset), family = poisson, 
      data = medpar)
  
  Coefficients:
               Estimate Std. Error  z value Pr(>|z|)    
  (Intercept) -3.151968   0.025880 -121.790   <2e-16 ***
  hmo         -0.001274   0.023735   -0.054    0.957    
  white       -0.289697   0.027284  -10.618   <2e-16 ***
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  (Dispersion parameter for poisson family taken to be 1)
  
      Null deviance: 11942  on 1494  degrees of freedom
  Residual deviance: 11837  on 1492  degrees of freedom
  AIC: 17558
  
  Number of Fisher Scoring iterations: 5
  
  > 
  > # How to obtain null deviance
  > 
  > irls(los ~ 1,
  +      family = "poisson",
  +      link = "log",
  +      offset = loset,
  +      data = medpar)$deviance
  [1] 11942.18
  > 
  > # cf
  > 
  > summary(glm.rpoi)$null.deviance
  [1] 11942.18
  > 
  > 
  > ## IDENTITY POISSON ----------------------------old-----------
  > 
  > irls.ipoi <- irls(los ~ hmo + white,
  +                   family = "poisson",
  +                   link = "identity",
  +                   data = medpar)
  > 
  > glm.ipoi <- glm(los ~ hmo + white,
  +                 family = poisson(link=identity),
  +                 data = medpar)
  > 
  > irls.ipoi
  
  Call:  irls(formula = los ~ hmo + white, data = medpar, family = "poisson", 
      link = "identity")
  
  Coefficients:
  X(Intercept)          Xhmo        Xwhite  
        11.926        -1.305        -2.036  
  
  Degrees of Freedom: 1494 Total (i.e. Null);  1492 Residual
  Null Deviance:	    NA 
  Residual Deviance: 8814 	AIC: 14530
  > glm.ipoi
  
  Call:  glm(formula = los ~ hmo + white, family = poisson(link = identity), 
      data = medpar)
  
  Coefficients:
  (Intercept)          hmo        white  
       11.926       -1.305       -2.036  
  
  Degrees of Freedom: 1494 Total (i.e. Null);  1492 Residual
  Null Deviance:	    8901 
  Residual Deviance: 8814 	AIC: 14530
  > 
  > summary(irls.ipoi)
  
  Call:
  irls(formula = los ~ hmo + white, data = medpar, family = "poisson", 
      link = "identity")
  
  Deviance Residuals:
     Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  -4.1104 -1.7492 -0.6218 -0.3074  0.9429 18.9468 
  
  Coefficients (all in linear predictor):
              Estimate    SE     Z        p   LCL    UCL
  (Intercept)    11.93 0.305 39.04        0 11.33 12.525
  hmo            -1.31 0.211 -6.19 5.84e-10 -1.72 -0.892
  white          -2.04 0.317 -6.43 1.24e-10 -2.66 -1.416
  
  Null deviance: NA  on  1494 d.f.
  Residual deviance: 8813.721  on  1492 d.f.
  Null Pearson:  on  1494 d.f.
  Residual Pearson: 11550.57  on  1492 d.f.
  Dispersion:
  AIC:  14534.87
  
  Number of optimizer iterations:  3 
  
  > summary(glm.ipoi)
  
  Call:
  glm(formula = los ~ hmo + white, family = poisson(link = identity), 
      data = medpar)
  
  Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
  (Intercept)  11.9263     0.3055  39.041  < 2e-16 ***
  hmo          -1.3053     0.2107  -6.195 5.84e-10 ***
  white        -2.0365     0.3165  -6.434 1.24e-10 ***
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  (Dispersion parameter for poisson family taken to be 1)
  
      Null deviance: 8901.1  on 1494  degrees of freedom
  Residual deviance: 8813.7  on 1492  degrees of freedom
  AIC: 14535
  
  Number of Fisher Scoring iterations: 4
  
  > 
  > ## NEGATIVE BINOMIAL (NB2) ------------------------------------------
  > 
  > irls.nb2 <- irls(los ~ hmo + white,
  +                  family = "negBinomial",
  +                  link = "log",
  +                  a = 0.5,
  +                  data = medpar)
  > 
  > glm.nb2 <- glm.nb(los ~ hmo + white,
  +                   data = medpar)
  > 
  > irls.nb2
  
  Call:  irls(formula = los ~ hmo + white, data = medpar, family = "negBinomial", 
      link = "log", a = 0.5)
  
  Coefficients:
  X(Intercept)          Xhmo        Xwhite  
        2.4810       -0.1405       -0.1897  
  
  Degrees of Freedom: 1494 Total (i.e. Null);  1492 Residual
  Null Deviance:	    NA 
  Residual Deviance: 1533 	AIC: 9705
  > glm.nb2
  
  Call:  glm.nb(formula = los ~ hmo + white, data = medpar, init.theta = 2.063546582, 
      link = log)
  
  Coefficients:
  (Intercept)          hmo        white  
       2.4810      -0.1405      -0.1897  
  
  Degrees of Freedom: 1494 Total (i.e. Null);  1492 Residual
  Null Deviance:	    1586 
  Residual Deviance: 1571 	AIC: 9706
  > 
  > summary(irls.nb2)
  
  Call:
  irls(formula = los ~ hmo + white, data = medpar, family = "negBinomial", 
      link = "log", a = 0.5)
  
  Deviance Residuals:
     Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  -2.0641 -0.8389 -0.2623 -0.2363  0.3724  5.4397 
  
  Coefficients (all in linear predictor):
              Estimate     SE     Z         p    LCL     UCL
  (Intercept)    2.481 0.0681 36.46 5.33e-291  2.348  2.6144
  hmo           -0.141 0.0553 -2.54    0.0111 -0.249 -0.0321
  white         -0.190 0.0711 -2.67   0.00766 -0.329 -0.0503
  
  Null deviance: NA  on  1494 d.f.
  Residual deviance: 1532.901  on  1492 d.f.
  Null Pearson:  on  1494 d.f.
  Residual Pearson: 1930.22  on  1492 d.f.
  Dispersion:
  AIC:  9704.601
  
  Number of optimizer iterations:  4 
  
  > summary(glm.nb2)
  
  Call:
  glm.nb(formula = los ~ hmo + white, data = medpar, init.theta = 2.063546582, 
      link = log)
  
  Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
  (Intercept)  2.48102    0.06715  36.947  < 2e-16 ***
  hmo         -0.14051    0.05463  -2.572  0.01011 *  
  white       -0.18971    0.07020  -2.702  0.00689 ** 
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  (Dispersion parameter for Negative Binomial(2.0635) family taken to be 1)
  
      Null deviance: 1585.7  on 1494  degrees of freedom
  Residual deviance: 1570.7  on 1492  degrees of freedom
  AIC: 9706.1
  
  Number of Fisher Scoring iterations: 1
  
  
                Theta:  2.0635 
            Std. Err.:  0.0893 
  
   2 x log-likelihood:  -9698.0790 
  > 
  > ## The IRLS deviance is wrong because of the conflict between JLL functions.
  > 
  > ## RATE NEGATIVE BINOMIAL (NB2)--------------------------------
  > 
  > irls.rnb2 <- irls(los ~ hmo + white,
  +                   family = "negBinomial",
  +                   link = "log",
  +                   a = 0.5,
  +                   offset=loset,
  +                   data = medpar)
  > 
  > glm.rnb2 <- glm.nb(los ~ hmo + white + offset(loset),
  +                    data = medpar)
  > 
  > irls.rnb2
  
  Call:  irls(formula = los ~ hmo + white, data = medpar, family = "negBinomial", 
      link = "log", offset = loset, a = 0.5)
  
  Coefficients:
  X(Intercept)          Xhmo        Xwhite  
      -2.97035       0.07026      -0.26306  
  
  Degrees of Freedom: 1494 Total (i.e. Null);  1492 Residual
  Null Deviance:	    NA 
  Residual Deviance: 2201 	AIC: 10370
  > glm.rnb2
  
  Call:  glm.nb(formula = los ~ hmo + white + offset(loset), data = medpar, 
      init.theta = 1.383408218, link = log)
  
  Coefficients:
  (Intercept)          hmo        white  
     -2.95831      0.07137     -0.25826  
  
  Degrees of Freedom: 1494 Total (i.e. Null);  1492 Residual
  Null Deviance:	    1644 
  Residual Deviance: 1633 	AIC: 10290
  > 
  > summary(irls.rnb2)
  
  Call:
  irls(formula = los ~ hmo + white, data = medpar, family = "negBinomial", 
      link = "log", offset = loset, a = 0.5)
  
  Deviance Residuals:
     Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  -2.4812 -1.1600 -0.4493 -0.2985  0.4008  5.3324 
  
  Coefficients (all in linear predictor):
              Estimate     SE      Z        p     LCL    UCL
  (Intercept)  -2.9704 0.0684 -43.43        0 -3.1044 -2.836
  hmo           0.0703 0.0553   1.27    0.204 -0.0382  0.179
  white        -0.2631 0.0715  -3.68 0.000233 -0.4032 -0.123
  
  Null deviance: NA  on  1494 d.f.
  Residual deviance: 2201.291  on  1492 d.f.
  Null Pearson:  on  1494 d.f.
  Residual Pearson: 2875.733  on  1492 d.f.
  Dispersion:
  AIC:  10372.99
  
  Number of optimizer iterations:  5 
  
  > summary(glm.rnb2)
  
  Call:
  glm.nb(formula = los ~ hmo + white + offset(loset), data = medpar, 
      init.theta = 1.383408218, link = log)
  
  Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
  (Intercept) -2.95831    0.08027 -36.852  < 2e-16 ***
  hmo          0.07137    0.06463   1.104  0.26943    
  white       -0.25826    0.08388  -3.079  0.00208 ** 
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  (Dispersion parameter for Negative Binomial(1.3834) family taken to be 1)
  
      Null deviance: 1644.0  on 1494  degrees of freedom
  Residual deviance: 1633.5  on 1492  degrees of freedom
  AIC: 10285
  
  Number of Fisher Scoring iterations: 1
  
  
                Theta:  1.3834 
            Std. Err.:  0.0542 
  
   2 x log-likelihood:  -10277.0180 
  > 
  > #oset <- rep(1:5, each=299, times=1)*100
  > #loset <- log(oset)
  > #glm.rnb2 <- glm(los ~ hmo + white + offset(loset),
  > #                family = negBinomial(2),
  > #                data = medpar)
  > #summary(glm.rnb2)
  > 
  > ## The IRLS deviance is wrong because of the conflict between JLL functions.
  > 
  > ###  CANONICAL NEGATIVE BINOMIAL (NBC)-------------------------
  > irls.nbc <- irls(los ~ hmo + white,
  +                  family = "negBinomial",
  +                  link = "negbin",
  +                  a = 0.5,
  +                  data = medpar)
  > irls.nbc
  
  Call:  irls(formula = los ~ hmo + white, data = medpar, family = "negBinomial", 
      link = "negbin", a = 0.5)
  
  Coefficients:
  X(Intercept)          Xhmo        Xwhite  
      -0.15428      -0.02499      -0.02996  
  
  Degrees of Freedom: 1494 Total (i.e. Null);  1492 Residual
  Null Deviance:	    NA 
  Residual Deviance: 1533 	AIC: 9705
  > summary(irls.nbc)
  
  Call:
  irls(formula = los ~ hmo + white, data = medpar, family = "negBinomial", 
      link = "negbin", a = 0.5)
  
  Deviance Residuals:
     Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  -2.0664 -0.8391 -0.2621 -0.2363  0.3727  5.4402 
  
  Coefficients (all in linear predictor):
              Estimate      SE      Z        p     LCL      UCL
  (Intercept)   -0.154 0.00984 -15.68 2.02e-55 -0.1736 -0.13500
  hmo           -0.025 0.01018  -2.45   0.0141 -0.0450 -0.00504
  white         -0.030 0.01047  -2.86  0.00421 -0.0505 -0.00944
  
  Null deviance: NA  on  1494 d.f.
  Residual deviance: 1532.809  on  1492 d.f.
  Null Pearson:  on  1494 d.f.
  Residual Pearson: 1929.641  on  1492 d.f.
  Dispersion:
  AIC:  9704.509
  
  Number of optimizer iterations:  4 
  
  > 
  > 
  > ###  IDENTITY NEGATIVE BINOMIAL --------------------------
  > irls.nbi <- irls(los ~ hmo + white,
  +                  family = "negBinomial",
  +                  link = "identity",
  +                  a = 0.5,
  +                  data = medpar)
  > irls.nbi
  
  Call:  irls(formula = los ~ hmo + white, data = medpar, family = "negBinomial", 
      link = "identity", a = 0.5)
  
  Coefficients:
  X(Intercept)          Xhmo        Xwhite  
        11.911        -1.295        -2.022  
  
  Degrees of Freedom: 1494 Total (i.e. Null);  1492 Residual
  Null Deviance:	    NA 
  Residual Deviance: 1533 	AIC: 9705
  > summary(irls.nbi)
  
  Call:
  irls(formula = los ~ hmo + white, data = medpar, family = "negBinomial", 
      link = "identity", a = 0.5)
  
  Deviance Residuals:
     Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  -2.0614 -0.8392 -0.2624 -0.2363  0.3723  5.4393 
  
  Coefficients (all in linear predictor):
              Estimate    SE     Z        p   LCL    UCL
  (Intercept)    11.91 0.801 14.87 4.81e-50 10.34 13.480
  hmo            -1.29 0.492 -2.63  0.00853 -2.26 -0.330
  white          -2.02 0.825 -2.45   0.0143 -3.64 -0.404
  
  Null deviance: NA  on  1494 d.f.
  Residual deviance: 1533.024  on  1492 d.f.
  Null Pearson:  on  1494 d.f.
  Residual Pearson: 1930.99  on  1492 d.f.
  Dispersion:
  AIC:  9704.724
  
  Number of optimizer iterations:  4 
  
  > 
  > 
  > 
  > ##  BERNOULLI LOGIT --------------------------------------------
  > 
  > irls.logit <- irls(died ~ hmo + white,
  +                    family = "binomial",
  +                    link = "logit",
  +                    data = medpar)
  > glm.logit <- glm(died ~ hmo + white,
  +                  family = "binomial",
  +                  data = medpar)
  > 
  > irls.logit
  
  Call:  irls(formula = died ~ hmo + white, data = medpar, family = "binomial", 
      link = "logit")
  
  Coefficients:
  X(Intercept)          Xhmo        Xwhite  
      -0.92619      -0.01225       0.30339  
  
  Degrees of Freedom: 1494 Total (i.e. Null);  1492 Residual
  Null Deviance:	    NA 
  Residual Deviance: 1921 	AIC: 1927
  > glm.logit
  
  Call:  glm(formula = died ~ hmo + white, family = "binomial", data = medpar)
  
  Coefficients:
  (Intercept)          hmo        white  
     -0.92619     -0.01225      0.30339  
  
  Degrees of Freedom: 1494 Total (i.e. Null);  1492 Residual
  Null Deviance:	    1923 
  Residual Deviance: 1921 	AIC: 1927
  > 
  > summary(irls.logit)
  
  Call:
  irls(formula = died ~ hmo + white, data = medpar, family = "binomial", 
      link = "logit")
  
  Deviance Residuals:
     Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  -0.9268 -0.9268 -0.9222 -0.1002  1.4507  1.5929 
  
  Coefficients (all in linear predictor):
              Estimate    SE       Z       p     LCL    UCL
  (Intercept)  -0.9262 0.197 -4.6922 2.7e-06 -1.3131 -0.539
  hmo          -0.0122 0.149 -0.0822   0.934 -0.3041  0.280
  white         0.3034 0.205  1.4787   0.139 -0.0988  0.706
  
  Null deviance: NA  on  1494 d.f.
  Residual deviance: 1920.602  on  1492 d.f.
  Null Pearson:  on  1494 d.f.
  Residual Pearson: 1495.004  on  1492 d.f.
  Dispersion:
  AIC:  1926.602
  
  Number of optimizer iterations:  3 
  
  > summary(glm.logit)
  
  Call:
  glm(formula = died ~ hmo + white, family = "binomial", data = medpar)
  
  Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
  (Intercept) -0.92619    0.19739  -4.692  2.7e-06 ***
  hmo         -0.01225    0.14893  -0.082    0.934    
  white        0.30339    0.20518   1.479    0.139    
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  (Dispersion parameter for binomial family taken to be 1)
  
      Null deviance: 1922.9  on 1494  degrees of freedom
  Residual deviance: 1920.6  on 1492  degrees of freedom
  AIC: 1926.6
  
  Number of Fisher Scoring iterations: 4
  
  > 
  > ## BINOMIAL LOGIT -------------------------------------------
  > 
  > denom <- rep(1:5, each=299, times=1)*100
  > 
  > irls.blogit <- irls(los ~ hmo + white,
  +                     family = "binomial",
  +                     link = "logit",
  +                     m = denom,
  +                     data = medpar)
  > 
  > not.los <- denom - medpar$los
  > glm.blogit <- glm(cbind(los, not.los) ~ hmo + white,
  +                   family = "binomial",
  +                   data = medpar)
  > 
  > glm.blogit
  
  Call:  glm(formula = cbind(los, not.los) ~ hmo + white, family = "binomial", 
      data = medpar)
  
  Coefficients:
  (Intercept)          hmo        white  
    -3.108260    -0.001317    -0.300870  
  
  Degrees of Freedom: 1494 Total (i.e. Null);  1492 Residual
  Null Deviance:	    12540 
  Residual Deviance: 12430 	AIC: 18090
  > irls.blogit
  
  Call:  irls(formula = los ~ hmo + white, data = medpar, family = "binomial", 
      link = "logit", m = denom)
  
  Coefficients:
  X(Intercept)          Xhmo        Xwhite  
     -3.108260     -0.001317     -0.300870  
  
  Degrees of Freedom: 1494 Total (i.e. Null);  1492 Residual
  Null Deviance:	    NA 
  Residual Deviance: 12430 	AIC: 18090
  > 
  > summary(irls.blogit)
  
  Call:
  irls(formula = los ~ hmo + white, data = medpar, family = "binomial", 
      link = "logit", m = denom)
  
  Deviance Residuals:
      Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
  -5.48311 -2.07549 -0.52043 -0.03487  1.64470 16.78605 
  
  Coefficients (all in linear predictor):
              Estimate     SE         Z        p     LCL    UCL
  (Intercept) -3.10826 0.0265 -117.5078        0 -3.1601 -3.056
  hmo         -0.00132 0.0241   -0.0546    0.956 -0.0486  0.046
  white       -0.30087 0.0279  -10.7954 3.62e-27 -0.3555 -0.246
  
  Null deviance: NA  on  1494 d.f.
  Residual deviance: 12433.1  on  1492 d.f.
  Null Pearson:  on  1494 d.f.
  Residual Pearson: 17223.7  on  1492 d.f.
  Dispersion:
  AIC:  18085.71
  
  Number of optimizer iterations:  5 
  
  > summary(glm.blogit)
  
  Call:
  glm(formula = cbind(los, not.los) ~ hmo + white, family = "binomial", 
      data = medpar)
  
  Coefficients:
               Estimate Std. Error  z value Pr(>|z|)    
  (Intercept) -3.108260   0.026452 -117.508   <2e-16 ***
  hmo         -0.001317   0.024132   -0.055    0.956    
  white       -0.300870   0.027870  -10.795   <2e-16 ***
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  (Dispersion parameter for binomial family taken to be 1)
  
      Null deviance: 12542  on 1494  degrees of freedom
  Residual deviance: 12433  on 1492  degrees of freedom
  AIC: 18086
  
  Number of Fisher Scoring iterations: 5
  
  > 
  > ## BERNOULLI PROBIT --------------------------------------------
  > 
  > irls.probit <- irls(died ~ hmo + white,
  +                     family = "binomial",
  +                     link = "probit",
  +                     data = medpar)
  > glm.probit <- glm(died ~ hmo + white,
  +                   family = binomial(link=probit),
  +                   data = medpar)
  > 
  > 
  > irls.probit
  
  Call:  irls(formula = died ~ hmo + white, data = medpar, family = "binomial", 
      link = "probit")
  
  Coefficients:
  X(Intercept)          Xhmo        Xwhite  
     -0.571875     -0.007304      0.184210  
  
  Degrees of Freedom: 1494 Total (i.e. Null);  1492 Residual
  Null Deviance:	    NA 
  Residual Deviance: 1921 	AIC: 1927
  > glm.probit
  
  Call:  glm(formula = died ~ hmo + white, family = binomial(link = probit), 
      data = medpar)
  
  Coefficients:
  (Intercept)          hmo        white  
    -0.571875    -0.007304     0.184210  
  
  Degrees of Freedom: 1494 Total (i.e. Null);  1492 Residual
  Null Deviance:	    1923 
  Residual Deviance: 1921 	AIC: 1927
  > 
  > summary(irls.probit)
  
  Call:
  irls(formula = died ~ hmo + white, data = medpar, family = "binomial", 
      link = "probit")
  
  Deviance Residuals:
     Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  -0.9268 -0.9268 -0.9223 -0.1002  1.4507  1.5928 
  
  Coefficients (all in linear predictor):
              Estimate     SE     Z        p     LCL    UCL
  (Intercept)  -0.5719 0.1184 -4.83 1.37e-06 -0.8040 -0.340
  hmo          -0.0073 0.0913 -0.08    0.936 -0.1862  0.172
  white         0.1842 0.1233  1.49    0.135 -0.0575  0.426
  
  Null deviance: NA  on  1494 d.f.
  Residual deviance: 1920.602  on  1492 d.f.
  Null Pearson:  on  1494 d.f.
  Residual Pearson: 1495.003  on  1492 d.f.
  Dispersion:
  AIC:  1926.602
  
  Number of optimizer iterations:  3 
  
  > summary(glm.probit)
  
  Call:
  glm(formula = died ~ hmo + white, family = binomial(link = probit), 
      data = medpar)
  
  Coefficients:
               Estimate Std. Error z value Pr(>|z|)    
  (Intercept) -0.571875   0.118412  -4.830 1.37e-06 ***
  hmo         -0.007304   0.091279  -0.080    0.936    
  white        0.184210   0.123302   1.494    0.135    
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  (Dispersion parameter for binomial family taken to be 1)
  
      Null deviance: 1922.9  on 1494  degrees of freedom
  Residual deviance: 1920.6  on 1492  degrees of freedom
  AIC: 1926.6
  
  Number of Fisher Scoring iterations: 4
  
  > 
  > ## BINOMIAL PROBIT------------------------------------------
  > 
  > denom <- rep(1:5, each=299, times=1)*100
  > 
  > irls.bprobit <- irls(los ~ hmo + white,
  +                      family = "binomial",
  +                      link = "probit",
  +                      m = denom,
  +                      data = medpar,
  +                      verbose = 1)
  1 -1.67458 0.007496065 -0.1476654 3158.497 
  2 -1.717767 1.202034e-05 -0.1334726 -82.23842 
  3 -1.719439 -0.0004795235 -0.1325953 -0.07489042 
  4 -1.719442 -0.0004810955 -0.1325929 -1.164735e-07 
  > 
  > not.los <- denom - medpar$los
  > glm.bprobit <- glm(cbind(los, not.los) ~ hmo + white,
  +                    family = binomial(link=probit),
  +                    data = medpar)
  > 
  > 
  > irls.bprobit
  
  Call:  irls(formula = los ~ hmo + white, data = medpar, family = "binomial", 
      link = "probit", m = denom, verbose = 1)
  
  Coefficients:
  X(Intercept)          Xhmo        Xwhite  
    -1.7194417    -0.0004811    -0.1325929  
  
  Degrees of Freedom: 1494 Total (i.e. Null);  1492 Residual
  Null Deviance:	    NA 
  Residual Deviance: 12430 	AIC: 18090
  > glm.bprobit
  
  Call:  glm(formula = cbind(los, not.los) ~ hmo + white, family = binomial(link = probit), 
      data = medpar)
  
  Coefficients:
  (Intercept)          hmo        white  
   -1.7194417   -0.0004811   -0.1325929  
  
  Degrees of Freedom: 1494 Total (i.e. Null);  1492 Residual
  Null Deviance:	    12540 
  Residual Deviance: 12430 	AIC: 18090
  > 
  > summary(irls.bprobit)
  
  Call:
  irls(formula = los ~ hmo + white, data = medpar, family = "binomial", 
      link = "probit", m = denom, verbose = 1)
  
  Deviance Residuals:
      Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
  -5.48304 -2.07541 -0.52031 -0.03487  1.64481 16.78622 
  
  Coefficients (all in linear predictor):
               Estimate     SE         Z        p     LCL    UCL
  (Intercept) -1.719442 0.0119 -144.4845        0 -1.7428 -1.696
  hmo         -0.000481 0.0104   -0.0461    0.963 -0.0209  0.020
  white       -0.132593 0.0125  -10.6167 2.49e-26 -0.1571 -0.108
  
  Null deviance: NA  on  1494 d.f.
  Residual deviance: 12433.1  on  1492 d.f.
  Null Pearson:  on  1494 d.f.
  Residual Pearson: 17223.7  on  1492 d.f.
  Dispersion:
  AIC:  18085.71
  
  Number of optimizer iterations:  4 
  
  > summary(glm.bprobit)
  
  Call:
  glm(formula = cbind(los, not.los) ~ hmo + white, family = binomial(link = probit), 
      data = medpar)
  
  Coefficients:
                Estimate Std. Error  z value Pr(>|z|)    
  (Intercept) -1.7194417  0.0119003 -144.487   <2e-16 ***
  hmo         -0.0004811  0.0104386   -0.046    0.963    
  white       -0.1325929  0.0124888  -10.617   <2e-16 ***
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  (Dispersion parameter for binomial family taken to be 1)
  
      Null deviance: 12542  on 1494  degrees of freedom
  Residual deviance: 12433  on 1492  degrees of freedom
  AIC: 18086
  
  Number of Fisher Scoring iterations: 4
  
  > 
  > ## BERNOULLI CLOGLOG --------------------------------------------
  > 
  > irls.cloglog <- irls(died ~ hmo + white,
  +                      family = "binomial",
  +                      link = "cloglog",
  +                      data = medpar,
  +                      verbose = 1)
  1 -1.082039 -0.01000365 0.2369577 -2.25514 
  2 -1.097596 -0.01045469 0.2524736 -0.008836906 
  3 -1.097676 -0.01048139 0.2525578 -2.899562e-07 
  > 
  > glm.cloglog <- glm(died ~ hmo + white,
  +                    family = binomial(link=cloglog),
  +                    data = medpar)
  > 
  > irls.cloglog
  
  Call:  irls(formula = died ~ hmo + white, data = medpar, family = "binomial", 
      link = "cloglog", verbose = 1)
  
  Coefficients:
  X(Intercept)          Xhmo        Xwhite  
      -1.09768      -0.01048       0.25256  
  
  Degrees of Freedom: 1494 Total (i.e. Null);  1492 Residual
  Null Deviance:	    NA 
  Residual Deviance: 1921 	AIC: 1927
  > glm.cloglog
  
  Call:  glm(formula = died ~ hmo + white, family = binomial(link = cloglog), 
      data = medpar)
  
  Coefficients:
  (Intercept)          hmo        white  
     -1.09768     -0.01048      0.25256  
  
  Degrees of Freedom: 1494 Total (i.e. Null);  1492 Residual
  Null Deviance:	    1923 
  Residual Deviance: 1921 	AIC: 1927
  > 
  > summary(irls.cloglog)
  
  Call:
  irls(formula = died ~ hmo + white, data = medpar, family = "binomial", 
      link = "cloglog", verbose = 1)
  
  Deviance Residuals:
     Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  -0.9268 -0.9268 -0.9220 -0.1002  1.4507  1.5929 
  
  Coefficients (all in linear predictor):
              Estimate    SE       Z        p     LCL    UCL
  (Intercept)  -1.0977 0.168 -6.5407 6.12e-11 -1.4266 -0.769
  hmo          -0.0105 0.122 -0.0862    0.931 -0.2487  0.228
  white         0.2526 0.174  1.4524    0.146 -0.0883  0.593
  
  Null deviance: NA  on  1494 d.f.
  Residual deviance: 1920.601  on  1492 d.f.
  Null Pearson:  on  1494 d.f.
  Residual Pearson: 1495.005  on  1492 d.f.
  Dispersion:
  AIC:  1926.601
  
  Number of optimizer iterations:  3 
  
  > summary(glm.cloglog)
  
  Call:
  glm(formula = died ~ hmo + white, family = binomial(link = cloglog), 
      data = medpar)
  
  Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
  (Intercept) -1.09768    0.16783  -6.540 6.13e-11 ***
  hmo         -0.01048    0.12154  -0.086    0.931    
  white        0.25256    0.17390   1.452    0.146    
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  (Dispersion parameter for binomial family taken to be 1)
  
      Null deviance: 1922.9  on 1494  degrees of freedom
  Residual deviance: 1920.6  on 1492  degrees of freedom
  AIC: 1926.6
  
  Number of Fisher Scoring iterations: 5
  
  > 
  > 
  > ## BINOMIAL CLOGLOG
  > 
  > denom <- rep(1:5, each=299, times=1)*100
  > 
  > irls.bcloglog <- irls(los ~ hmo + white,
  +                       family = "binomial",
  +                       link = "cloglog",
  +                       m = denom,
  +                       data = medpar,
  +                       verbose = 1)
  1 -2.966468 0.04118462 -0.3198788 3417.894 
  2 -3.118026 0.004789742 -0.2984934 -340.1226 
  3 -3.13012 -0.001252375 -0.295279 -1.588474 
  4 -3.130191 -0.001323089 -0.2952463 -4.923369e-05 
  5 -3.130191 -0.001323088 -0.2952463 0 
  > 
  > not.los <- denom - medpar$los
  > glm.bcloglog <- glm(cbind(los, not.los) ~ hmo + white,
  +                     family = binomial(link=cloglog),
  +                     data = medpar)
  > 
  > irls.bcloglog
  
  Call:  irls(formula = los ~ hmo + white, data = medpar, family = "binomial", 
      link = "cloglog", m = denom, verbose = 1)
  
  Coefficients:
  X(Intercept)          Xhmo        Xwhite  
     -3.130191     -0.001323     -0.295246  
  
  Degrees of Freedom: 1494 Total (i.e. Null);  1492 Residual
  Null Deviance:	    NA 
  Residual Deviance: 12430 	AIC: 18090
  > glm.bcloglog
  
  Call:  glm(formula = cbind(los, not.los) ~ hmo + white, family = binomial(link = cloglog), 
      data = medpar)
  
  Coefficients:
  (Intercept)          hmo        white  
    -3.130191    -0.001323    -0.295246  
  
  Degrees of Freedom: 1494 Total (i.e. Null);  1492 Residual
  Null Deviance:	    12540 
  Residual Deviance: 12430 	AIC: 18090
  > 
  > summary(irls.bcloglog)
  
  Call:
  irls(formula = los ~ hmo + white, data = medpar, family = "binomial", 
      link = "cloglog", m = denom, verbose = 1)
  
  Deviance Residuals:
      Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
  -5.48312 -2.07550 -0.52044 -0.03486  1.64468 16.78602 
  
  Coefficients (all in linear predictor):
              Estimate     SE         Z        p     LCL     UCL
  (Intercept) -3.13019 0.0259 -120.9384        0 -3.1809 -3.0795
  hmo         -0.00132 0.0237   -0.0557    0.956 -0.0478  0.0452
  white       -0.29525 0.0273  -10.8205 2.75e-27 -0.3487 -0.2418
  
  Null deviance: NA  on  1494 d.f.
  Residual deviance: 12433.1  on  1492 d.f.
  Null Pearson:  on  1494 d.f.
  Residual Pearson: 17223.69  on  1492 d.f.
  Dispersion:
  AIC:  18085.71
  
  Number of optimizer iterations:  5 
  
  > summary(glm.bcloglog)
  
  Call:
  glm(formula = cbind(los, not.los) ~ hmo + white, family = binomial(link = cloglog), 
      data = medpar)
  
  Coefficients:
               Estimate Std. Error  z value Pr(>|z|)    
  (Intercept) -3.130191   0.025883 -120.938   <2e-16 ***
  hmo         -0.001323   0.023736   -0.056    0.956    
  white       -0.295246   0.027286  -10.821   <2e-16 ***
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  (Dispersion parameter for binomial family taken to be 1)
  
      Null deviance: 12542  on 1494  degrees of freedom
  Residual deviance: 12433  on 1492  degrees of freedom
  AIC: 18086
  
  Number of Fisher Scoring iterations: 5
  
  > 
  > 
  > ## GAMMA CANONICAL--------------------------------------
  > #irls.gam <- irls(los ~ hmo + white,
  > #                 family = "gamma",
  > #                 link = "inverse",
  > #                 data = medpar)
  > 
  > #glm.gam <- glm(los ~ hmo + white,
  > #               family = Gamma,
  > #               data = medpar)
  > 
  > #irls.gam
  > #glm.gam
  > 
  > #summary(irls.gam)
  > #summary(glm.gam)
  > 
  > 
  > ## INVERSE GAUSSIAN CANONICAL--------------------------------------
  > #irls.ivg <- irls(los ~ hmo + white,
  > #                 family = "inv_gauss",
  > #                 link = "inverse2",
  > #                 data = medpar)
  > #glm.ivg <- glm(los ~ hmo + white,
  > #               family = inverse.gaussian,
  > #               data = medpar)
  > 
  > #irls.ivg
  > #glm.ivg
  > 
  > #summary(irls.ivg)
  > #summary(glm.ivg)
  > 
  > 
  > proc.time()
     user  system elapsed 
     0.79    0.29    1.03 
Running the tests in 'tests/test_ml_g.R' failed.
Complete output:
  > 
  > 
  > library(msme)
  Loading required package: MASS
  Loading required package: lattice
  > 
  > data(ufc)
  > ufc <- na.omit(ufc)
  > 
  > ufc.g.reg <- ml_g(height.m ~ dbh.cm, data = ufc)
  > 
  > ufc.g.lm <- lm(height.m ~ dbh.cm, data = ufc)
  > 
  > ufc.g.reg
  
  Call:
  ml_g(formula = height.m ~ dbh.cm, data = ufc)
  
  No coefficients
  
  > ufc.g.lm
  
  Call:
  lm(formula = height.m ~ dbh.cm, data = ufc)
  
  Coefficients:
  (Intercept)       dbh.cm  
      12.6757       0.3126  
  
  > 
  > summary(ufc.g.reg)
  
  Call:
  ml_g(formula = height.m ~ dbh.cm, data = ufc)
  
  Residuals:
      Min      1Q  Median      3Q     Max 
  -33.525  -2.863   0.132   2.851  13.320 
  
  Coefficients:
              Estimate       SE        Z      LCL    UCL
  (Intercept) 12.67708  0.56271 22.52853 11.57416 13.780
  dbh.cm       0.31256  0.01384 22.57707  0.28543  0.340
  Log Sigma    1.59521  0.03577 44.60167  1.52511  1.665
  
  Residual standard error: 4.929 on 391 degrees of freedom
  
  > summary(ufc.g.lm)
  
  Call:
  lm(formula = height.m ~ dbh.cm, data = ufc)
  
  Residuals:
      Min      1Q  Median      3Q     Max 
  -33.526  -2.862   0.132   2.851  13.321 
  
  Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
  (Intercept) 12.67570    0.56406   22.47   <2e-16 ***
  dbh.cm       0.31259    0.01388   22.52   <2e-16 ***
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  Residual standard error: 4.941 on 389 degrees of freedom
  Multiple R-squared:  0.566,	Adjusted R-squared:  0.5649 
  F-statistic: 507.4 on 1 and 389 DF,  p-value: < 2.2e-16
  
  > 
  > coef(ufc.g.reg)
  NULL
  > 
  > logLik(ufc.g.reg)
  'log Lik.' NaN (df=)
  > logLik(ufc.g.lm)
  'log Lik.' -1178.469 (df=3)
  > 
  > AIC(ufc.g.reg)
  numeric(0)
  > AIC(ufc.g.lm)
  [1] 2362.938
  > 
  > 
  > 
  > 
  > 
  > proc.time()
     user  system elapsed 
     0.29    0.07    0.34 
Running the tests in 'tests/test_ml_glm2.R' failed.
Complete output:
  > library(msme)
  Loading required package: MASS
  Loading required package: lattice
  > 
  > # library(msme, lib.loc="lib")
  > 
  > library(MASS)
  > 
  > data(medpar)
  > data(ufc)
  > 
  > denom <- rep(1:5, each=299, times=1)*100   # m : binomial denominator w medpar
  > oset <- rep(1:5, each=299, times=1)*100    # offset Poisson, NB, offset w medpar
  > loset <- log(oset)                         # log of oset
  > 
  > ## Normal
  > 
  > ufc <- na.omit(ufc)
  > ml.g <- ml_glm2(height.m ~ dbh.cm,
  +                 formula2 = ~1,
  +                 data = ufc,
  +                 family = "normal",
  +                 mean.link = "identity",
  +                 scale.link = "log_s")
  > 
  > lm.g <- lm(height.m ~ dbh.cm,
  +                 data = ufc)
  >                 
  > ml.g
  
  Call:  ml_glm2(formula1 = height.m ~ dbh.cm, formula2 = ~1, data = ufc, 
      family = "normal", mean.link = "identity", scale.link = "log_s")
  
  Coefficients:
    (Intercept)         dbh.cm  (Intercept)_s  
        12.6757         0.3126         1.5950  
  
  Degrees of Freedom: 389 Total (i.e. Null);  388 Residual
  Null Deviance:	    21890 
  Residual Deviance: 9498 	AIC: 2363
  > lm.g
  
  Call:
  lm(formula = height.m ~ dbh.cm, data = ufc)
  
  Coefficients:
  (Intercept)       dbh.cm  
      12.6757       0.3126  
  
  > 
  > summary(ml.g)
  
  Call:
  ml_glm2(formula1 = height.m ~ dbh.cm, formula2 = ~1, data = ufc, 
      family = "normal", mean.link = "identity", scale.link = "log_s")
  
  Deviance Residuals:
     Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  -33.526  -2.862   0.132   0.000   2.851  13.321 
  
  Pearson Residuals:
       Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
  -2.072894 -0.260070  0.009491 -0.006867  0.232057  1.074036 
  
  Coefficients (all in linear predictor):
                Estimate     SE    Z         p    LCL   UCL
  (Intercept)     12.676 0.5626 22.5 2.12e-112 11.573 13.78
  dbh.cm           0.313 0.0138 22.6 6.37e-113  0.285  0.34
  (Intercept)_s    1.595 0.0358 44.6         0  1.525  1.67
  
  Null deviance: 21885.77  on  389 d.f.
  Residual deviance: 9497.697  on  388 d.f.
  Null Pearson: 154.5773  on  389 d.f.
  Residual Pearson: 66.69231  on  388 d.f.
  Dispersion: 0.1718874
  AIC:  2362.938
  
  Number of optimizer iterations:  47 
  
  > summary(lm.g)
  
  Call:
  lm(formula = height.m ~ dbh.cm, data = ufc)
  
  Residuals:
      Min      1Q  Median      3Q     Max 
  -33.526  -2.862   0.132   2.851  13.321 
  
  Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
  (Intercept) 12.67570    0.56406   22.47   <2e-16 ***
  dbh.cm       0.31259    0.01388   22.52   <2e-16 ***
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  Residual standard error: 4.941 on 389 degrees of freedom
  Multiple R-squared:  0.566,	Adjusted R-squared:  0.5649 
  F-statistic: 507.4 on 1 and 389 DF,  p-value: < 2.2e-16
  
  > 
  > ## NEGATIVE BINOMIAL (NB2) ------------------------------------------
  > 
  > glm.nb2 <- glm.nb(los ~ hmo + white,
  +                   data = medpar)
  > 
  > ml.nb2 <- ml_glm2(los ~ hmo + white,
  +                     formula2 = ~1,
  +                     data = medpar,
  +                     family = "negBinomial",
  +                     mean.link = "log",
  +                     scale.link = "log_s")
  > 
  > glm.nb2
  
  Call:  glm.nb(formula = los ~ hmo + white, data = medpar, init.theta = 2.063546582, 
      link = log)
  
  Coefficients:
  (Intercept)          hmo        white  
       2.4810      -0.1405      -0.1897  
  
  Degrees of Freedom: 1494 Total (i.e. Null);  1492 Residual
  Null Deviance:	    1586 
  Residual Deviance: 1571 	AIC: 9706
  > ml.nb2
  
  Call:  ml_glm2(formula1 = los ~ hmo + white, formula2 = ~1, data = medpar, 
      family = "negBinomial", mean.link = "log", scale.link = "log_s")
  
  Coefficients:
    (Intercept)            hmo          white  (Intercept)_s  
         2.4810        -0.1405        -0.1897        -0.7244  
  
  Degrees of Freedom: 1493 Total (i.e. Null);  1491 Residual
  Null Deviance:	    1586 
  Residual Deviance: 1571 	AIC: 9706
  > 
  > summary(glm.nb2)
  
  Call:
  glm.nb(formula = los ~ hmo + white, data = medpar, init.theta = 2.063546582, 
      link = log)
  
  Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
  (Intercept)  2.48102    0.06715  36.947  < 2e-16 ***
  hmo         -0.14051    0.05463  -2.572  0.01011 *  
  white       -0.18971    0.07020  -2.702  0.00689 ** 
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  (Dispersion parameter for Negative Binomial(2.0635) family taken to be 1)
  
      Null deviance: 1585.7  on 1494  degrees of freedom
  Residual deviance: 1570.7  on 1492  degrees of freedom
  AIC: 9706.1
  
  Number of Fisher Scoring iterations: 1
  
  
                Theta:  2.0635 
            Std. Err.:  0.0893 
  
   2 x log-likelihood:  -9698.0790 
  > summary(ml.nb2)
  
  Call:
  ml_glm2(formula1 = los ~ hmo + white, formula2 = ~1, data = medpar, 
      family = "negBinomial", mean.link = "log", scale.link = "log_s")
  
  Deviance Residuals:
     Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  -2.0870 -0.8490 -0.2657 -0.2385  0.3774  5.5163 
  
  Pearson Residuals:
       Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
  -0.625347 -0.361972 -0.129771  0.000002  0.213921  7.293970 
  
  Coefficients (all in linear predictor):
                Estimate     SE      Z         p    LCL     UCL
  (Intercept)      2.481 0.0671  36.97 3.69e-299  2.349  2.6126
  hmo             -0.141 0.0547  -2.57    0.0101 -0.248 -0.0334
  white           -0.190 0.0702  -2.70    0.0069 -0.327 -0.0521
  (Intercept)_s   -0.724 0.0433 -16.74  6.73e-63 -0.809 -0.6396
  
  Null deviance: 1585.658  on  1493 d.f.
  Residual deviance: 1570.678  on  1491 d.f.
  Null Pearson: 554.4386  on  1493 d.f.
  Residual Pearson: 535.8635  on  1491 d.f.
  Dispersion: 0.3593987
  AIC:  9706.079
  
  Number of optimizer iterations:  56 
  
  > 
  > ## RATE NEGATIVE BINOMIAL (NB2)--------------------------------
  > 
  > glm.rnb2 <- glm.nb(los ~ hmo + white + offset(loset),
  +                    data = medpar)
  > 
  > ml.rnb2 <- ml_glm2(los ~ hmo + white,
  +                     formula2 = ~1,
  +                     data = medpar,
  +                   offset=loset,
  +                     family = "negBinomial",
  +                     mean.link = "log",
  +                     scale.link = "inverse_s")
  There were 35 warnings (use warnings() to see them)
  > 
  > ml.rnb2
  
  Call:  ml_glm2(formula1 = los ~ hmo + white, formula2 = ~1, data = medpar, 
      family = "negBinomial", mean.link = "log", scale.link = "inverse_s", 
      offset = loset)
  
  Coefficients:
    (Intercept)            hmo          white  (Intercept)_s  
       -2.95831        0.07137       -0.25826        1.38341  
  
  Degrees of Freedom: 1493 Total (i.e. Null);  1491 Residual
  Null Deviance:	    1644 
  Residual Deviance: 1633 	AIC: 10290
  > glm.rnb2
  
  Call:  glm.nb(formula = los ~ hmo + white + offset(loset), data = medpar, 
      init.theta = 1.383408218, link = log)
  
  Coefficients:
  (Intercept)          hmo        white  
     -2.95831      0.07137     -0.25826  
  
  Degrees of Freedom: 1494 Total (i.e. Null);  1492 Residual
  Null Deviance:	    1644 
  Residual Deviance: 1633 	AIC: 10290
  > 
  > summary(glm.rnb2)
  
  Call:
  glm.nb(formula = los ~ hmo + white + offset(loset), data = medpar, 
      init.theta = 1.383408218, link = log)
  
  Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
  (Intercept) -2.95831    0.08027 -36.852  < 2e-16 ***
  hmo          0.07137    0.06463   1.104  0.26943    
  white       -0.25826    0.08388  -3.079  0.00208 ** 
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  (Dispersion parameter for Negative Binomial(1.3834) family taken to be 1)
  
      Null deviance: 1644.0  on 1494  degrees of freedom
  Residual deviance: 1633.5  on 1492  degrees of freedom
  AIC: 10285
  
  Number of Fisher Scoring iterations: 1
  
  
                Theta:  1.3834 
            Std. Err.:  0.0542 
  
   2 x log-likelihood:  -10277.0180 
  > summary(ml.rnb2)
  
  Call:
  ml_glm2(formula1 = los ~ hmo + white, formula2 = ~1, data = medpar, 
      family = "negBinomial", mean.link = "log", scale.link = "inverse_s", 
      offset = loset)
  
  Deviance Residuals:
     Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  -2.1562 -1.0031 -0.3981 -0.2739  0.3196  4.5472 
  
  Pearson Residuals:
      Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
  -0.79744 -0.54361 -0.26121  0.02789  0.25739  8.97718 
  
  Coefficients (all in linear predictor):
                Estimate     SE      Z         p     LCL     UCL
  (Intercept)    -2.9583 0.0814 -36.34 4.34e-289 -3.1179 -2.7987
  hmo             0.0714 0.0660   1.08     0.279 -0.0579  0.2007
  white          -0.2583 0.0851  -3.03   0.00242 -0.4251 -0.0914
  (Intercept)_s   1.3834 0.0544  25.45 7.12e-143  1.2769  1.4900
  
  Null deviance: 1643.981  on  1493 d.f.
  Residual deviance: 1633.494  on  1491 d.f.
  Null Pearson: 1217.899  on  1493 d.f.
  Residual Pearson: 1212.726  on  1491 d.f.
  Dispersion: 0.8133639
  AIC:  10285.02
  
  Number of optimizer iterations:  32 
  
  > 
  > 
  > 
  > ## GAMMA CANONICAL--------------------------------------
  > #irls.gam <- irls(los ~ hmo + white,
  > #                 family = "gamma",
  > #                 link = "inverse",
  > #                 data = medpar)
  > 
  > #glm.gam <- glm(los ~ hmo + white,
  > #               family = Gamma,
  > #               data = medpar)
  > 
  > #ml.gam <- ml_glm2(los ~ hmo + white,
  > #                    formula2 = ~1, 
  > #                    data = medpar,
  > #                    family = "gamma",
  > #                    mean.link = "inverse",
  > #                    scale.link = "inverse_s")
  > 
  > #irls.gam
  > #glm.gam
  > #ml.gam
  > 
  > #summary(irls.gam)
  > #summary(glm.gam)
  > #summary(ml.gam)
  > 
  > 
  > ## INVERSE GAUSSIAN CANONICAL--------------------------------------
  > #irls.ivg <- irls(los ~ hmo + white,
  > #                 family = "inv_gauss",
  > #                 link = "inverse2",
  > #                 data = medpar)
  > #glm.ivg <- glm(los ~ hmo + white,
  > #               family = inverse.gaussian,
  > #               data = medpar)
  > 
  > #irls.ivg
  > #glm.ivg
  > 
  > #summary(irls.ivg)
  > #summary(glm.ivg)
  > 
  > 
  > proc.time()
     user  system elapsed 
     1.56    0.18    1.73 
Running the tests in 'tests/test_nbinomial.R' failed.
Complete output:
  > # nbinomial.r    msme package
  > #                Joseph M Hilbe and Andrew P Robinson
  > #                Methods of Statistical Model Estimation
  > #                Chapman & Hall/CRC  2013
  > # Examples of use for the nbinomial.r function for maximum likelihood 
  > # NB2 and heterogeneous NB2 regression. The default "family" is 
  > # "nb2", which estimates parameters and the dispersion parameter using 
  > # a direct relationship between the dispersion and Poisson variance. 
  > # The default NB variance is mu + alpha*mu^2. Use the "negBinomial" family
  > # to produce the inverted dispersion, which is given by R's glm.nb function
  > 
  > library(msme)
  Loading required package: MASS
  Loading required package: lattice
  > 
  > # library(msme, lib.loc="lib")
  > 
  > data(medpar)
  > 
  > # TRADITIONAL NB REGRESSION WITH ALPHA
  > mynb1 <- nbinomial(los ~ hmo + white, data=medpar)
  There were 34 warnings (use warnings() to see them)
  > summary(mynb1)
  
  Call:
  ml_glm2(formula1 = formula1, formula2 = formula2, data = data, 
      family = family, mean.link = mean.link, scale.link = scale.link, 
      offset = offset, start = start, verbose = verbose)
  
  Deviance Residuals:
     Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  -2.0870 -0.8490 -0.2657 -0.2385  0.3774  5.5163 
  
  Pearson Residuals:
       Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
  -1.215580 -0.689385 -0.249473  0.000002  0.411242 14.021982 
  
  Coefficients (all in linear predictor):
                Estimate     SE     Z         p    LCL     UCL
  (Intercept)      2.481 0.0671 36.97 3.69e-299  2.349  2.6126
  hmo             -0.141 0.0547 -2.57    0.0101 -0.248 -0.0334
  white           -0.190 0.0702 -2.70    0.0069 -0.327 -0.0521
  (Intercept)_s    0.485 0.0210 23.11 3.87e-118  0.443  0.5257
  
  Null deviance: 1585.656  on  1493 d.f.
  Residual deviance: 1570.675  on  1491 d.f.
  Null Pearson: 2048.131  on  1493 d.f.
  Residual Pearson: 1980.982  on  1491 d.f.
  Dispersion: 1.328627
  AIC:  9706.079
  
  Number of optimizer iterations:  68 
  
  > 
  > # DISPERSION STATISTIC A
  > mynb1$dispersion
  [1] 1.328627
  > 
  > # INCIDENCE RATE RATIOS
  > exp(mynb1$coefficients)
    (Intercept)           hmo         white (Intercept)_s 
     11.9534384     0.8689159     0.8272014     1.6235311 
  > 
  > # IRR SEs USING DELTA METHOD
  > exp(mynb1$coefficients)*mynb1$se.beta.hat
    (Intercept)           hmo         white (Intercept)_s 
     0.80221877    0.04748915    0.05809000    0.03404792 
  > 
  > #IRR CONFIDENCE INTERVALS
  > mynb1$coefficients - 1.96*mynb1$se.beta.hat
    (Intercept)           hmo         white (Intercept)_s 
      2.3494795    -0.2476294    -0.3273476     0.4434993 
  > mynb1$coefficients + 1.96*mynb1$se.beta.hat
    (Intercept)           hmo         white (Intercept)_s 
     2.61255843   -0.03338840   -0.05206662    0.52570768 
  > 
  > # TRADITIONAL NB -- SUMMARY INCLUDED IN FUNCTION CALL
  > summary(mynb1_5 <- nbinomial(los ~ hmo + white, data=medpar))
  
  Call:
  ml_glm2(formula1 = formula1, formula2 = formula2, data = data, 
      family = family, mean.link = mean.link, scale.link = scale.link, 
      offset = offset, start = start, verbose = verbose)
  
  Deviance Residuals:
     Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  -2.0870 -0.8490 -0.2657 -0.2385  0.3774  5.5163 
  
  Pearson Residuals:
       Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
  -1.215580 -0.689385 -0.249473  0.000002  0.411242 14.021982 
  
  Coefficients (all in linear predictor):
                Estimate     SE     Z         p    LCL     UCL
  (Intercept)      2.481 0.0671 36.97 3.69e-299  2.349  2.6126
  hmo             -0.141 0.0547 -2.57    0.0101 -0.248 -0.0334
  white           -0.190 0.0702 -2.70    0.0069 -0.327 -0.0521
  (Intercept)_s    0.485 0.0210 23.11 3.87e-118  0.443  0.5257
  
  Null deviance: 1585.656  on  1493 d.f.
  Residual deviance: 1570.675  on  1491 d.f.
  Null Pearson: 2048.131  on  1493 d.f.
  Residual Pearson: 1980.982  on  1491 d.f.
  Dispersion: 1.328627
  AIC:  9706.079
  
  Number of optimizer iterations:  68 
  
  There were 34 warnings (use warnings() to see them)
  > 
  > # TRADITIONAL NB -- SHOWING ALL OPTIONS
  > mynb2 <- nbinomial(los ~ hmo + white,
  +                     formula2 = ~ 1,
  +                     data = medpar,
  +                     family = "nb2",
  +                     mean.link = "log",
  +                     scale.link = "inverse_s")
  There were 34 warnings (use warnings() to see them)
  > summary(mynb2)
  
  Call:
  ml_glm2(formula1 = formula1, formula2 = formula2, data = data, 
      family = family, mean.link = mean.link, scale.link = scale.link, 
      offset = offset, start = start, verbose = verbose)
  
  Deviance Residuals:
     Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  -2.0870 -0.8490 -0.2657 -0.2385  0.3774  5.5163 
  
  Pearson Residuals:
       Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
  -1.215580 -0.689385 -0.249473  0.000002  0.411242 14.021982 
  
  Coefficients (all in linear predictor):
                Estimate     SE     Z         p    LCL     UCL
  (Intercept)      2.481 0.0671 36.97 3.69e-299  2.349  2.6126
  hmo             -0.141 0.0547 -2.57    0.0101 -0.248 -0.0334
  white           -0.190 0.0702 -2.70    0.0069 -0.327 -0.0521
  (Intercept)_s    0.485 0.0210 23.11 3.87e-118  0.443  0.5257
  
  Null deviance: 1585.656  on  1493 d.f.
  Residual deviance: 1570.675  on  1491 d.f.
  Null Pearson: 2048.131  on  1493 d.f.
  Residual Pearson: 1980.982  on  1491 d.f.
  Dispersion: 1.328627
  AIC:  9706.079
  
  Number of optimizer iterations:  68 
  
  > 
  > # R GLM.NB - LIKE INVERTED DISPERSION BASED M
  > mynb3 <- nbinomial(los ~ hmo + white,
  +                     formula2 = ~ 1,
  +                     data = medpar,
  +                     family = "negBinomial",
  +                     mean.link = "log",
  +                     scale.link = "inverse_s")
  There were 15 warnings (use warnings() to see them)
  > summary(mynb3)
  
  Call:
  ml_glm2(formula1 = formula1, formula2 = formula2, data = data, 
      family = family, mean.link = mean.link, scale.link = scale.link, 
      offset = offset, start = start, verbose = verbose)
  
  Deviance Residuals:
     Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  -2.0870 -0.8490 -0.2657 -0.2385  0.3774  5.5163 
  
  Pearson Residuals:
       Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
  -0.625347 -0.361972 -0.129771  0.000002  0.213921  7.293970 
  
  Coefficients (all in linear predictor):
                Estimate     SE     Z         p    LCL     UCL
  (Intercept)      2.481 0.0671 36.97 3.69e-299  2.349  2.6126
  hmo             -0.141 0.0547 -2.57    0.0101 -0.248 -0.0334
  white           -0.190 0.0702 -2.70    0.0069 -0.327 -0.0521
  (Intercept)_s    2.064 0.0893 23.11 3.89e-118  1.889  2.2386
  
  Null deviance: 1585.659  on  1493 d.f.
  Residual deviance: 1570.678  on  1491 d.f.
  Null Pearson: 554.4386  on  1493 d.f.
  Residual Pearson: 535.8634  on  1491 d.f.
  Dispersion: 0.3593987
  AIC:  9706.079
  
  Number of optimizer iterations:  64 
  
  > 
  > # R GLM.NB-TYPE INVERTED DISPERSON --THETA ; WITH DEFAULTS
  > mynb4 <- nbinomial(los ~ hmo + white, family="negBinomial", data =medpar)
  There were 15 warnings (use warnings() to see them)
  > summary(mynb4)
  
  Call:
  ml_glm2(formula1 = formula1, formula2 = formula2, data = data, 
      family = family, mean.link = mean.link, scale.link = scale.link, 
      offset = offset, start = start, verbose = verbose)
  
  Deviance Residuals:
     Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  -2.0870 -0.8490 -0.2657 -0.2385  0.3774  5.5163 
  
  Pearson Residuals:
       Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
  -0.625347 -0.361972 -0.129771  0.000002  0.213921  7.293970 
  
  Coefficients (all in linear predictor):
                Estimate     SE     Z         p    LCL     UCL
  (Intercept)      2.481 0.0671 36.97 3.69e-299  2.349  2.6126
  hmo             -0.141 0.0547 -2.57    0.0101 -0.248 -0.0334
  white           -0.190 0.0702 -2.70    0.0069 -0.327 -0.0521
  (Intercept)_s    2.064 0.0893 23.11 3.89e-118  1.889  2.2386
  
  Null deviance: 1585.659  on  1493 d.f.
  Residual deviance: 1570.678  on  1491 d.f.
  Null Pearson: 554.4386  on  1493 d.f.
  Residual Pearson: 535.8634  on  1491 d.f.
  Dispersion: 0.3593987
  AIC:  9706.079
  
  Number of optimizer iterations:  64 
  
  > 
  > # HETEROGENEOUS NB; DISPERSION PARAMETERIZED
  > mynb5 <- nbinomial(los ~ hmo + white,
  +                     formula2 = ~ hmo + white,
  +                     data = medpar,
  +                     family = "negBinomial",
  +                     mean.link = "log",
  +                     scale.link = "log_s")
  > summary(mynb5)
  
  Call:
  ml_glm2(formula1 = formula1, formula2 = formula2, data = data, 
      family = family, mean.link = mean.link, scale.link = scale.link, 
      offset = offset, start = start, verbose = verbose)
  
  Deviance Residuals:
     Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  -2.0798 -0.9203 -0.2630 -0.2383  0.3729  5.4506 
  
  Pearson Residuals:
       Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
  -0.652743 -0.340323 -0.131515  0.000018  0.216499  7.385571 
  
  Coefficients (all in linear predictor):
                Estimate     SE      Z         p    LCL      UCL
  (Intercept)     2.4788 0.0691 35.885 5.18e-282  2.343  2.61414
  hmo            -0.1403 0.0511 -2.746   0.00602 -0.240 -0.04018
  white          -0.1873 0.0720 -2.599   0.00934 -0.328 -0.04607
  (Intercept)_s  -0.6345 0.1397 -4.542  5.56e-06 -0.908 -0.36072
  hmo_s          -0.2432 0.1279 -1.902    0.0572 -0.494  0.00746
  white_s        -0.0634 0.1471 -0.431     0.667 -0.352  0.22490
  
  Null deviance: 1477.833  on  1493 d.f.
  Residual deviance: 1572.623  on  1489 d.f.
  Null Pearson: 603.9427  on  1493 d.f.
  Residual Pearson: 540.9986  on  1489 d.f.
  Dispersion: 0.3633301
  AIC:  9706.119
  
  Number of optimizer iterations:  74 
  
  > 
  > # SAVED STATISTICS FROM NBINOMIAL2.R
  > 
  > 
  > # USE OF "nbh" for glm.nb-like dipersion
  > # nbinomial now has "negBinomial" as the family for glm.nb
  > 
  > 
  > # TRADITIONAL NB2 REGRESSION. USING DEFAULT OPTIONS
  > tnb1 <- nbinomial(los ~ hmo + white, data=medpar)
  There were 34 warnings (use warnings() to see them)
  > summary(tnb1)
  
  Call:
  ml_glm2(formula1 = formula1, formula2 = formula2, data = data, 
      family = family, mean.link = mean.link, scale.link = scale.link, 
      offset = offset, start = start, verbose = verbose)
  
  Deviance Residuals:
     Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  -2.0870 -0.8490 -0.2657 -0.2385  0.3774  5.5163 
  
  Pearson Residuals:
       Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
  -1.215580 -0.689385 -0.249473  0.000002  0.411242 14.021982 
  
  Coefficients (all in linear predictor):
                Estimate     SE     Z         p    LCL     UCL
  (Intercept)      2.481 0.0671 36.97 3.69e-299  2.349  2.6126
  hmo             -0.141 0.0547 -2.57    0.0101 -0.248 -0.0334
  white           -0.190 0.0702 -2.70    0.0069 -0.327 -0.0521
  (Intercept)_s    0.485 0.0210 23.11 3.87e-118  0.443  0.5257
  
  Null deviance: 1585.656  on  1493 d.f.
  Residual deviance: 1570.675  on  1491 d.f.
  Null Pearson: 2048.131  on  1493 d.f.
  Residual Pearson: 1980.982  on  1491 d.f.
  Dispersion: 1.328627
  AIC:  9706.079
  
  Number of optimizer iterations:  68 
  
  > 
  > 
  > # TRADITIONAL NB2 REGRESSION. PROVIDING EXPLICIT OPTIONS
  > tnb2 <- nbinomial(los ~ hmo + white,
  +                     formula2 = ~ 1,
  +                     data = medpar,
  +                     family = "nb2",
  +                     mean.link = "log",
  +                     scale.link = "inverse_s")
  There were 34 warnings (use warnings() to see them)
  > 
  > summary(tnb2)
  
  Call:
  ml_glm2(formula1 = formula1, formula2 = formula2, data = data, 
      family = family, mean.link = mean.link, scale.link = scale.link, 
      offset = offset, start = start, verbose = verbose)
  
  Deviance Residuals:
     Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  -2.0870 -0.8490 -0.2657 -0.2385  0.3774  5.5163 
  
  Pearson Residuals:
       Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
  -1.215580 -0.689385 -0.249473  0.000002  0.411242 14.021982 
  
  Coefficients (all in linear predictor):
                Estimate     SE     Z         p    LCL     UCL
  (Intercept)      2.481 0.0671 36.97 3.69e-299  2.349  2.6126
  hmo             -0.141 0.0547 -2.57    0.0101 -0.248 -0.0334
  white           -0.190 0.0702 -2.70    0.0069 -0.327 -0.0521
  (Intercept)_s    0.485 0.0210 23.11 3.87e-118  0.443  0.5257
  
  Null deviance: 1585.656  on  1493 d.f.
  Residual deviance: 1570.675  on  1491 d.f.
  Null Pearson: 2048.131  on  1493 d.f.
  Residual Pearson: 1980.982  on  1491 d.f.
  Dispersion: 1.328627
  AIC:  9706.079
  
  Number of optimizer iterations:  68 
  
  > 
  > # PEARSON CHI2 STATISTIC
  > tnb1$pearson
  [1] 1980.982
  > 
  > 
  > 
  > 
  > proc.time()
     user  system elapsed 
     5.92    0.39    6.31