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) > > 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 ------------------------------------------------- > > glm.poi <- glm(los ~ hmo + white, + family = "poisson", + data = medpar) > > ml.poi <- ml_glm(los ~ hmo + white, + family = "poisson", + link = "log", + data = medpar) Warning message: In y - y.hat : Recycling array of length 1 in vector-array arithmetic is deprecated. Use c() or as.vector() instead. > > ml.poi Call: ml_glm(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: 8901 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(ml.poi) Call: ml_glm(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: 8901.134 on 1494 d.f. Residual deviance: 8812.942 on 1492 d.f. AIC: 14534.09 Number of optimizer iterations: 54 > 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 > > ## RATE POISSON > > glm.rpoi <- glm(los ~ hmo + white + offset(loset), + family = poisson, + data = medpar) > ml.rpoi <- ml_glm(los ~ hmo + white, + family = "poisson", link = "log", + offset = loset, + data = medpar) There were 50 or more warnings (use warnings() to see the first 50) > > ml.rpoi Call: ml_glm(formula = los ~ hmo + white, data = medpar, family = "poisson", link = "log", offset = loset) Coefficients: X(Intercept) Xhmo Xwhite -3.151968 -0.001274 -0.289697 Degrees of Freedom: 1494 Total (i.e. Null); 1492 Residual Null Deviance: 11940 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(ml.rpoi) Call: ml_glm(formula = los ~ hmo + white, data = medpar, family = "poisson", link = "log", offset = loset) Deviance Residuals: Min. 1st Qu. Median Mean 3rd Qu. Max. -5.41201 -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: 11942.18 on 1494 d.f. Residual deviance: 11837.32 on 1492 d.f. AIC: 17558.47 Number of optimizer iterations: 51 > 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 > > ## IDENTITY POISSON ----------------------------old----------- > > glm.ipoi <- glm(los ~ hmo + white, + family = poisson(link=identity), + data = medpar) > > ml.ipoi <- ml_glm(los ~ hmo + white, + family = "poisson", link = "identity", + data = medpar) Warning message: In y - y.hat : Recycling array of length 1 in vector-array arithmetic is deprecated. Use c() or as.vector() instead. > > 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 > ml.ipoi Call: ml_glm(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: 8901 Residual Deviance: 8814 AIC: 14530 > > 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 > summary(ml.ipoi) Call: ml_glm(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.306 39.04 0 11.33 12.525 hmo -1.31 0.211 -6.19 5.86e-10 -1.72 -0.892 white -2.04 0.317 -6.43 1.31e-10 -2.66 -1.415 Null deviance: 8901.134 on 1494 d.f. Residual deviance: 8813.721 on 1492 d.f. AIC: 14534.87 Number of optimizer iterations: 15 > > ## BERNOULLI LOGIT -------------------------------------------- > > glm.logit <- glm(died ~ hmo + white, + family = "binomial", + data = medpar) > ml.logit <- ml_glm(died ~ hmo + white, + data = medpar, + family = "bernoulli", + link = "logit1") Warning message: In y - y.hat : Recycling array of length 1 in vector-array arithmetic is deprecated. Use c() or as.vector() instead. > > ml.logit Call: ml_glm(formula = died ~ hmo + white, data = medpar, family = "bernoulli", link = "logit1") Coefficients: X(Intercept) Xhmo Xwhite -0.92619 -0.01225 0.30339 Degrees of Freedom: 1494 Total (i.e. Null); 1492 Residual Null Deviance: 1923 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(ml.logit) Call: ml_glm(formula = died ~ hmo + white, data = medpar, family = "bernoulli", link = "logit1") 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.4786 0.139 -0.0988 0.706 Null deviance: 1922.865 on 1494 d.f. Residual deviance: 1920.602 on 1492 d.f. AIC: 1926.602 Number of optimizer iterations: 32 > 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 > > ## BERNOULLI PROBIT -------------------------------------------- > > ml.probit <- ml_glm(died ~ hmo + white, + family = "bernoulli", + link = "probit1", + data = medpar) Warning message: In y - y.hat : Recycling array of length 1 in vector-array arithmetic is deprecated. Use c() or as.vector() instead. > glm.probit <- glm(died ~ hmo + white, + family = binomial(link=probit), + data = medpar) > > > ml.probit Call: ml_glm(formula = died ~ hmo + white, data = medpar, family = "bernoulli", link = "probit1") Coefficients: X(Intercept) Xhmo Xwhite -0.571875 -0.007304 0.184210 Degrees of Freedom: 1494 Total (i.e. Null); 1492 Residual Null Deviance: 1923 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(ml.probit) Call: ml_glm(formula = died ~ hmo + white, data = medpar, family = "bernoulli", link = "probit1") 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.0574 0.426 Null deviance: 1922.865 on 1494 d.f. Residual deviance: 1920.602 on 1492 d.f. AIC: 1926.602 Number of optimizer iterations: 44 > 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 > > ## BERNOULLI CLOGLOG -------------------------------------------- > > ml.cloglog <- ml_glm(died ~ hmo + white, + family = "bernoulli", + link = "cloglog1", + data = medpar, + verbose = 1) $par X(Intercept) Xhmo Xwhite -1.09767633 -0.01048153 0.25255795 $value [1] -960.3007 $counts function gradient 50 14 $convergence [1] 0 $message NULL $hessian X(Intercept) Xhmo Xwhite X(Intercept) -505.42738 -80.80018 -469.75872 Xhmo -80.80018 -80.80018 -77.57093 Xwhite -469.75872 -77.57093 -469.75872 Warning message: In y - y.hat : Recycling array of length 1 in vector-array arithmetic is deprecated. Use c() or as.vector() instead. > > glm.cloglog <- glm(died ~ hmo + white, + family = binomial(link=cloglog), + data = medpar) > > ml.cloglog Call: ml_glm(formula = died ~ hmo + white, data = medpar, family = "bernoulli", link = "cloglog1", 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: 1923 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(ml.cloglog) Call: ml_glm(formula = died ~ hmo + white, data = medpar, family = "bernoulli", link = "cloglog1", 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.5416 6.09e-11 -1.4266 -0.769 hmo -0.0105 0.122 -0.0862 0.931 -0.2487 0.228 white 0.2526 0.174 1.4522 0.146 -0.0883 0.593 Null deviance: 1922.865 on 1494 d.f. Residual deviance: 1920.601 on 1492 d.f. AIC: 1926.601 Number of optimizer iterations: 50 > 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 > > > proc.time() user system elapsed 1.04 0.21 1.26