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. > 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.6253468 -0.3619718 -0.1297711 0.0000022 0.2139206 7.2939702 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: 53 > > ## 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.54360 -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: 35 > > > > ## 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.60 0.25 1.84