[6562] in s-news-athena
daemon@ATHENA.MIT.EDU (YONG CHANG )
Fri Mar 3 16:36:52 1995
From: ychang@welchlink.welch.jhu.edu (YONG CHANG )
To: s-news@utstat.toronto.edu
Date: Fri, 3 Mar 1995 16:07:30 -0500 (EST)
Due to an error in my rebuilt alias list, the following message hadn't
reached S newsgroup as expected. Now I'm trying to send it over again,
and I apologize that it'd taken so long for me to realize what it'd gone
wrong.
------------------------------------------------------------------------------
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
}
