[R-sig-eco] offset in Logistic regression

[Andrew Muir]

Dear list,

I am attempting to include an offset variable (Julian Day) in a series of 
logistic models using the glm() function in R.
I want to account for Juliann day of sampling in each of these models 
because sample collection was temporally stratified across the entire year.
The models work fine without the offset variable, but when I include it I 
get the following error messages.

> #Full model (Global hypothesis)
> M1= glm(Presence~Lat + Depth + Cliff + Temp_var + offset(J_Day), 
> family=binomial(link=logit))
Warning messages:
1: In glm.fit(x = X, y = Y, weights = weights, start = start, etastart = 
etastart,  :
  algorithm did not converge
2: In glm.fit(x = X, y = Y, weights = weights, start = start, etastart = 
etastart,  :
  fitted probabilities numerically 0 or 1 occurred
3: In glm.fit(x = X[, "(Intercept)", drop = FALSE], y = Y, weights = 
weights,  :
  fitted probabilities numerically 0 or 1 occurred
> summary (M1)

Warning messages:
1: In method(x = x[, varseq <= i, drop = FALSE], y = object$y, weights = 
object$prior.weights,  :
  algorithm did not converge
2: In method(x = x[, varseq <= i, drop = FALSE], y = object$y, weights = 
object$prior.weights,  :
  fitted probabilities numerically 0 or 1 occurred
3: In method(x = x[, varseq <= i, drop = FALSE], y = object$y, weights = 
object$prior.weights,  :
  fitted probabilities numerically 0 or 1 occurred
4: In method(x = x[, varseq <= i, drop = FALSE], y = object$y, weights = 
object$prior.weights,  :
  fitted probabilities numerically 0 or 1 occurred
>


I spent considerable time on the internet trying to figure out what the 
problem is and I came across the following bug report from 2002:

[Rd] glm start/offset bugs (PR#1422)

presnell at stat.ufl.edu presnell at stat.ufl.edu
Fri, 29 Mar 2002 04:04:03 +0100 (MET)

  a.. Previous message: [Rd] pos argument to library
  b.. Next message: [Rd] glm start/offset bugs (PR#1422)
  c.. Messages sorted by: [ date ] [ thread ] [ subject ] [ author ]

--------------------------------------------------------------------------------

--fupGvOGOQM

Content-Type: text/plain; charset=us-ascii

Content-Description: message body and .signature

Content-Transfer-Encoding: 7bit





There's a simple bug in the handling of the start and offset arguments

in glm and glm.fit.  The bug exists in the latest development version

of R (version information below), but it appears that glm.R has not

been touched much lately, so the bug affects at least the most recent

stable release of R.



Here is a simple example:



> data(ships, package=MASS)

> ships.glm <- glm(incidents ~ type + year + period + offset(log(service)),

+                family=poisson, data=ships, subset=service != 0)

> glm(incidents ~ type + year + period + offset(log(service)),

+     family=poisson, data=ships, subset=service != 0, 
start=coef(ships.glm))



or more simply



> update(ships.glm, start=coef(ships.glm))



The problem is caused by a bad initialization of etastart in glm.fit()

when an offset is present.  Fixing this and retrying the example above

then reveals another bug, this time in glm(), that is tickled only

when non-NULL values are given to both the offset and start arguments.



I've attached to the bottom of this message a simple patch to

src/library/base/R/glm.R that I believe correctly repairs these bugs.



Essentially the same change is needed in Berwin Turlach's Wed, 27 Feb

2002 version of glm.fit() reported on the R bug tracking site as

"Models/1331".  (The change to glm() is still required as well of

course.)  As an aside, I'm wondering if there has been any thought

given to adopting Berwin's patches.  Though I must admit that I

haven't taken the time to check carefully through his changes, I know

Berwin to be a very careful and skillful programmer, and I suspect

that he is the only person to have looked very closely at those

particular details of glm.fit in recent months (years?).  Also, some

of the simpler changes that he made seem to be needed just to cleanup

the code.  Would it be useful if a number of us began using his

version as an informal test?  After discovering this bug earlier

today, I have already begun to do so and so far have not encountered

any problems.




I sent an e-mail to Dr. Presnell (the person that wrote the bug report), but 
to-date I have not received a response from him.
So, my questions to the list are:
1. does anyone know if this bug has been fixed in R 2.6.1? (Dr. Presnell 
included a patch for this problem in his bug report, but I don't know if I 
should try to install the patch)
2. perhaps my problem is not with the glm.fit function, but a problem with 
the structure of my data or with my use of the offset command, if so, can 
anyone provide suggestions?

Thanks in advance,

AMM.

Andrew M. Muir
Purdue University,
Department of Forestry and Natural Resources,
Pfendler Hall of Agriculture, Rm G004
715 W. State St.
West Lafayette, IN  47907
Voice: 765-494-9597
FAX: 765-496-2422
email: ammuir at purdue.edu