Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

glm starting values not accepted log-link

Tags:

r

offset

glm

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) ?

like image 275
Andi Avatar asked Nov 21 '11 12:11

Andi


2 Answers

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.

like image 102
Ben Bolker Avatar answered Nov 08 '22 16:11

Ben Bolker


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
like image 27
James Avatar answered Nov 08 '22 15:11

James