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 ------------------------------------------------- > > 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.52042 -0.03486 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.81 0.31 1.10