I want to run a Gaussian GLM with a log link and an offset. The following problems arise:
y <- c(1,1,0,0)
t <- c(5,3,2,4)
No problem:
exp(coef(glm(y~1 + offset(log(t)), family=poisson)))
with family=gaussian
, starting values need to be specified, it works here:
exp(coef(glm(y~1, family=gaussian(link=log), start=0)))
but does not work here:
exp(coef(glm(y~1 + offset(log(t)), family=gaussian(link=log), start=0)))
Error in eval(expr, envir, enclos) : cannot find valid starting values: please specify some"
Does anyone see what's wrong (hopefully just in my coding) ?
Here are the results of some archaeology that explains what's going on, deep within the glm
function:
Debugging (with debug("glm")
) and stepping through the function shows that it fails at the following call:
if (length(offset) && attr(mt, "intercept") > 0L) {
fit$null.deviance <- eval(call(if (is.function(method)) "method" else method,
x = X[, "(Intercept)", drop = FALSE], y = Y, weights = weights,
offset = offset, family = family, control = control,
intercept = TRUE))$deviance
}
This is an attempt to calculate the null deviance for the model. It's only evaluated if there's an intercept term and an offset term (I'm not sure why; it may be that the default null deviance calculated by the previous call to glm
is wrong in that case and must be recalculated?). It calls glm.fit
(the default value of method
), but without starting values because these are usually unnecessary for the intercept-only model.
Now debugging inside glm.fit
to see what happens: we get (within a call to the family function, gaussian()
) to:
if (is.null(etastart) && is.null(start) && is.null(mustart) &&
((family$link == "inverse" && any(y == 0)) || (family$link ==
"log" && any(y <= 0))))
stop("cannot find valid starting values: please specify some")
and we see that because the starting values were not passed through, because a log link is used, and because some y
values are equal to zero, the fit fails. So this is a case that should happen if (and only if?) an offset and an intercept are both specified, a log link is used, and there are zero values in the response.
If you dump("glm",file="glmtemp.R")
; add the line
start = start[1], etastart = etastart[1], mustart = mustart[1],
to the call that fits the null deviance (i.e. the one shown above); and source("glmtemp.R")
, it seems to work OK ... I think this should be a reasonable general solution. If anyone wants to bring this issue up on the R development list, feel free.
I looks like start
isn't being recognised when offset
is present. You are trying to take the log of 0 in the y
values which is -Inf
. glm
obviously cannot deal with this when looking for a solution without being given some help by start
. A small perturbation in your y
values will permit a solution.
exp(coef(glm(I(y+.Machine$double.eps)~1 + offset(log(t)), family=gaussian(link=log))))
(Intercept)
0.1481481
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With