[R] Logit / ms

[Prof Brian Ripley]

On Fri, 22 Feb 2002, John Janmaat wrote:

> Hello,
>
> I am looking for a routine to do a logistic regression.  In the book
> "Modern Applied Statistics with S-PLUS" a function is described which
> uses the 'ms' command.  Is there an analogue for R, or an alternative
> approach that can accomplish the same thing?

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.

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