[6561] in s-news-athena
Summary on Poisson Regression
daemon@ATHENA.MIT.EDU (YONG CHANG )
Fri Mar 3 15:01:20 1995
From: ychang@welchlink.welch.jhu.edu (YONG CHANG )
To: s-news@utstat.toronto.edu (S Newsgroup)
Date: Fri, 3 Mar 1995 14:10:26 -0500 (EST)
Due to some system error, the following summary apparently had never
reached S newsgroup as expected. So I'm trying to post it again, and
hopefully it can help others running poisson regression in the old
version of Splus 3.0.
-------------------------------------------------------------------------------
This message is a followup on the message I posted two days ago regarding
some "unbelievable" null deviance I had encountered running Poisson regression
in Splus.
Prof. Brian Ripley and Trevor Hastie have offered help to me in a very timely
manner, and both pointed out that the error could be due to omission of offset
in the Null model in the old version 3.0 of Splus, which turned out to be
the case. I was using 3.0 version of Splus, instead of the latest 3.2(3.3?).
The following is the corrected script of "glm" function given by Prof. Ripley.
Here's the current (corrected) glm():
glm
function(formula = formula(data), family = gaussian, data = sys.parent(),
weights, subset, na.action, start = eta, control = glm.control(...),
method = "glm.fit", model = F, x = F, y = T, contrasts = NULL, ...)
{
call <- match.call()
m <- match.call(expand = F)
m$family <- m$method <- m$model <- m$x <- m$y <- m$control <- m$
contrasts <- m$... <- NULL
m[[1]] <- as.name("model.frame")
m <- eval(m, sys.parent())
Terms <- attr(m, "terms")
if(method == "model.frame")
return(m)
a <- attributes(m)
Y <- model.extract(m, response)
X <- model.matrix(Terms, m, contrasts)
w <- model.extract(m, weights)
if(!length(w))
w <- rep(1, nrow(m))
else if(any(w < 0))
stop("negative weights not allowed")
start <- model.extract(m, start)
offset <- model.extract(m, offset)
family <- as.family(family)
if(missing(method))
method <- attr(family, "method")
if(!is.null(method)) {
if(!exists(method, mode = "function"))
stop(paste("unimplemented method:", method))
}
else method <- "glm.fit"
glm.fitter <- get(method)
fit <- glm.fitter(x = X, y = Y, w = w, start = start, offset = offset,
family = family, maxit = control$maxit, epsilon = control$
# If an offset and intercept is present, iterations are needed to
# compute the Null deviance; these are done here, unless the model
# is NULL, in which case the computations have been done already
#
if(any(offset) && attr(Terms, "intercept")) {
null.deviance <- if(length(Terms)) glm.fitter(X[, "(Intercept)",
drop = F], Y, w, offset = offset, family =
family, maxit = control$maxit, epsilon =
control$epsilon, null.dev = NULL)$deviance
else fit$deviance
fit$null.deviance <- null.deviance
}
attr(fit, "class") <- c("glm", "lm")
fit$terms <- Terms
fit$formula <- as.vector(attr(Terms, "formula"))
fit$call <- call
if(model)
fit$model <- m
if(x)
fit$x <- X
if(!y)
fit$y <- NULL
fit
}
