[R] Logit / ms

[Prof Brian D Ripley]

Sorry, my mailer has lost all my answers ....

On Sun, 24 Feb 2002, Prof Brian D Ripley wrote:

> On Sat, 23 Feb 2002, Paul Johnson wrote:
>
> > Thanks for posting this. it is highly instructive!
> >
> > Can I ask follow ups? I ran this example after getting the bwt data as
> > illustrated in the example for birthwt in MASS.  It runs fine and gives
> > me the parameter estimates.

Let me remind people of the context: this is an illustration in the
optimization section of MASS the book, and we were asked for an R version.
It is *not* a function from the MASS library.  In particular, it is a
bare-bones function taking minimal space in the text, illustrating
optimization.  There are important comments on p.269 on the context.

All page references are to the third edition.

> > Question 1. the estimates are a little different from the glm estimates
> > obtained. The differences result from a change in optimization routines?
> >   Are these small differences typical?
> >
> > Here are the logitreg() numbers:
> >
> > (Intercept)         age         lwt   raceblack   raceother   smokeTRUE
> >   0.82304295 -0.03723343 -0.01565330  1.19240547  0.74067565  0.75551956
> >      ptdTRUE      htTRUE      uiTRUE        ftv1       ftv2+
> >   1.34374814  1.91317620  0.68020276 -0.43636831  0.17901477
> >
> >  > glm(low ~ . ,binomial, bwt)
> >
> > Call:  glm(formula = low ~ ., family = binomial, data = bwt)
> >
> > Coefficients:
> > (Intercept)          age          lwt    raceblack    raceother    smokeTRUE
> >      0.82271     -0.03722     -0.01565      1.19223      0.74051
> > 0.75537
> >      ptdTRUE       htTRUE       uiTRUE         ftv1        ftv2+
> >      1.34365      1.91297      0.68016     -0.43633      0.17894

Yes. glm has too sloppy a convergence criterion.  Turn epsilon down in
glm.control.  See p.216.

> > Question 2. Then I wondered "how do I do significance tests on those
> > estimates"?  In the glm results, I use summary(). But what of this
> > logitreg? I figure just to use t tests based on the asymptotic normality
> > of the b's, so I need standard errors.  To get them, it appears to me I
> > go into the logitreg function, and for optim I insert Hessian=TRUE, and
> > then I can torture the Hessian to get standard errors.

No, you *need* to do likelihood-ratio tests as t-tests are unreliable.
See p.225.

> > Question 3. when logitreg prints its output, the only diagnostic
> > information it gives is:
> > Residual Deviance: 195.4755
> >
> > I'm wondering what the user is supposed to conclude from that. Isn't it
> > the same as -2LL?  What benchmark do you use to say it is high or low?
> > In the olden days of graduate school, they ignore that, and instead look
> > for -2LLR to test that all the b's are jointly 0.

Same issues as for glm, so look `deviance' up in the book's index.

It is very unusual to fit a model expecting all the coefficients (even the
intercept!) to be zero.  That is totally implausible in the example, where
it corresponds to 50-50 prior odds of a baby having low birth weight, for
every baby (and `low' is defined below as the 10percentile in most SGA
studies).  I won't ask which graduate school misled you!


-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272860 (secr)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._