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