[R] Logit / ms

[Paul Johnson]

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.

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


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.

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.

pj


Prof Brian Ripley wrote:
> 
> Well, you want the upcoming fourth edition, which just happens to have a
> version for R:
> 
> logitreg <- function(x, y, wt = rep(1, length(y)),
>                intercept = T, start = rep(0, p), ...)
> {
>   fmin <- function(beta, X, y, w) {
>       p <- plogis(X %*% beta)
>       -sum(2 * w * ifelse(y, log(p), log(1-p)))
>   }
>   gmin <- function(beta, X, y, w) {
>       eta <- X %*% beta; p <- plogis(eta)
>       -2 * matrix(w *dlogis(eta) * ifelse(y, 1/p, -1/(1-p)), 1) %*% X
>   }
>   if(is.null(dim(x))) dim(x) <- c(length(x), 1)
>   dn <- dimnames(x)[[2]]
>   if(!length(dn)) dn <- paste("Var", 1:ncol(x), sep="")
>   p <- ncol(x) + intercept
>   if(intercept) {x <- cbind(1, x); dn <- c("(Intercept)", dn)}
>   if(is.factor(y)) y <- (unclass(y) != 1)
>   fit <- optim(start, fmin, gmin, X = x, y = y, w = wt,
>                method = "BFGS", ...)
>   names(fit$par) <- dn
>   cat("\nCoefficients:\n"); print(fit$par)
>   # R: use fit$value and fit$convergence
>   cat("\nResidual Deviance:", format(fit$value), "\n")
>   if(fit$convergence > 0)
>       cat("\nConvergence code:", fit$convergence, "\n")
>   invisible(fit)
> }
> 
> options(contrasts = c("contr.treatment", "contr.poly"))
> X <- model.matrix(terms(low ~ ., data=bwt), data = bwt)[, -1]
> logitreg(X, bwt$low)
> 
> 
> Don't totally overlook glm, though.
> 
> 



-- 
Paul E. Johnson                       email: pauljohn at ukans.edu
Dept. of Political Science            http://lark.cc.ku.edu/~pauljohn
University of Kansas                  Office: (785) 864-9086
Lawrence, Kansas 66045                FAX: (785) 864-5700

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