R Under development (unstable) (2025-04-26 r88181 ucrt) -- "Unsuffered Consequences" Copyright (C) 2025 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > # 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.2155800 -0.6893847 -0.2494731 0.0000021 0.4112425 14.0219821 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: 66 > > # 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.2155800 -0.6893847 -0.2494731 0.0000021 0.4112425 14.0219821 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: 66 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.2155800 -0.6893847 -0.2494731 0.0000021 0.4112425 14.0219821 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: 66 > > # 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.6253467 -0.3619717 -0.1297710 0.0000022 0.2139207 7.2939701 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.8635 on 1491 d.f. Dispersion: 0.3593987 AIC: 9706.079 Number of optimizer iterations: 58 > > # 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.6253467 -0.3619717 -0.1297710 0.0000022 0.2139207 7.2939701 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.8635 on 1491 d.f. Dispersion: 0.3593987 AIC: 9706.079 Number of optimizer iterations: 58 > > # 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.6527433 -0.3403226 -0.1315146 0.0000185 0.2164986 7.3855713 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.2155800 -0.6893847 -0.2494731 0.0000021 0.4112425 14.0219821 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: 66 > > > # 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.2155800 -0.6893847 -0.2494731 0.0000021 0.4112425 14.0219821 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: 66 > > # PEARSON CHI2 STATISTIC > tnb1$pearson [1] 1980.982 > > > > > proc.time() user system elapsed 5.01 0.60 5.60