I've noticed that glm and lrm give extremely different results if you
attempt to fit a saturated model to a dataset with zero cells. Consider,
for instance the data from, Agresti's Death Penalty example [0].
The crosstab table is:
, , PENALTY = NO
VIC
DEF BLACK WHITE
BLACK 97 52
WHITE 9 132
, , PENALTY = YES
VIC
DEF BLACK WHITE
BLACK 6 11
WHITE 0 19
Regression with an unsaturated model produces essentially
the same fit parameters with both glm and lrm. However,
if we try to fit a saturated model....
FITTING WITH GLM:> summary(glm(PENALTY~DEF*VIC,binomial))
Call:
glm(formula = PENALTY ~ DEF * VIC, family = binomial)
Deviance Residuals:
Min 1Q Median 3Q Max
-0.6195 -0.5186 -0.5186 -0.3465 2.3845
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.7830 0.4207 -6.615 3.71e-11 ***
DEFWHITE -4.7823 8.8981 -0.537 0.5910
VICWHITE 1.2296 0.5358 2.295 0.0217 *
DEFWHITE:VICWHITE 4.3973 8.9076 0.494 0.6216
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` '
1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 226.51 on 325 degrees of freedom
Residual deviance: 218.39 on 322 degrees of freedom
AIC: 226.39
Number of Fisher Scoring iterations: 6
FITTING WITH LRM:> lrm(PENALTY~DEF*VIC)
Logistic Regression Model
lrm(formula = PENALTY ~ DEF * VIC)
Frequencies of Responses
NO YES
290 36
Obs Max Deriv Model L.R. d.f. P C Dxy
326 0.002 8.13 3 0.0435 0.624 0.248
Gamma Tau-a R2 Brier
0.383 0.049 0.049 0.096
Coef S.E. Wald Z P
Intercept -2.783 0.4207 -6.62 0.0000
DEF=WHITE -5.490 20.8691 -0.26 0.7925
VIC=WHITE 1.230 0.5358 2.29 0.0217
DEF=WHITE * VIC=WHITE 5.105 20.8732 0.24 0.8068
If we fill in the remaining table cell with a dummy value, [1]
however, then glm and lrm produce essentially the same result.
Here's the glm result.
> summary(glm(PENALTY~DEF*VIC,binomial))
Call:
glm(formula = PENALTY ~ DEF * VIC, family = binomial)
Deviance Residuals:
Min 1Q Median 3Q Max
-0.6195 -0.5186 -0.5186 -0.3465 2.3845
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.7829 0.4192 -6.639 3.15e-11 ***
DEFWHITE 0.5857 1.1343 0.516 0.6056
VICWHITE 1.2296 0.5346 2.300 0.0215 *
DEFWHITE:VICWHITE -0.9707 1.2070 -0.804 0.4213
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` '
1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 230.90 on 326 degrees of freedom
Residual deviance: 224.88 on 323 degrees of freedom
AIC: 232.88
Number of Fisher Scoring iterations: 4
So, my question here is: is this normal behavior? If it
is, perhaps someone could speculate on why the results are
different.
Thanks,
-Ekr
[0] Agresti, A., "Categorical Data Analysis", Wiley 1990.
The data set can be found at http://www.rtfm.com/death.txt
[1] http://www.rtfm.com/death-filled-in.txt
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at
stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Eric Rescorla <ekr at rtfm.com> writes:> I've noticed that glm and lrm give extremely different results if you > attempt to fit a saturated model to a dataset with zero cells. Consider,...> Coefficients: > Estimate Std. Error z value Pr(>|z|) > (Intercept) -2.7830 0.4207 -6.615 3.71e-11 *** > DEFWHITE -4.7823 8.8981 -0.537 0.5910 > VICWHITE 1.2296 0.5358 2.295 0.0217 * > DEFWHITE:VICWHITE 4.3973 8.9076 0.494 0.6216 > --- > Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 > > (Dispersion parameter for binomial family taken to be 1) > > Null deviance: 226.51 on 325 degrees of freedom > Residual deviance: 218.39 on 322 degrees of freedom...> Obs Max Deriv Model L.R. d.f. P C Dxy > 326 0.002 8.13 3 0.0435 0.624 0.248 > Gamma Tau-a R2 Brier > 0.383 0.049 0.049 0.096 > > Coef S.E. Wald Z P > Intercept -2.783 0.4207 -6.62 0.0000 > DEF=WHITE -5.490 20.8691 -0.26 0.7925 > VIC=WHITE 1.230 0.5358 2.29 0.0217 > DEF=WHITE * VIC=WHITE 5.105 20.8732 0.24 0.8068These are *not* extremely different! In fact, they are essentially equivalent. The maximum likelihood estimates of the 2nd and 4th coefficients are theoretically infinite, and it is only a matter of when the two routines decide to stop the iterations. The 1st and 3rd coefficients are the same, so are their SEs, and also the sum of the 2nd and 4th is -.39 in both cases. The likelihood ratio in the first model is 226.51-218.39 = 8.12 and in the second it is given as 8.13. [And where did lrm() come from? It's not in the standard packages, and I'm not going to wade through 150+ CRAN packages to locate it...] -- O__ ---- Peter Dalgaard Blegdamsvej 3 c/ /'_ --- Dept. of Biostatistics 2200 Cph. N (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk) FAX: (+45) 35327907 -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
On Thu, 24 Oct 2002, Eric Rescorla wrote:> I've noticed that glm and lrm give extremely different results if you > attempt to fit a saturated model to a dataset with zero cells. Consider, > for instance the data from, Agresti's Death Penalty example [0].The MLEs are infinite, so there are bound to be problems. Both methods will stop short of infinity, when the coefficients or likelihood stop changing, but where they stop will depend on the exact tolerance and stopping rule. They aren't actually that different in terms of fitted values. For example, the zero cell has a fitted probability of 0.00025 and 0.0005 in the two models. -thomas -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._