Curious finding in MASS:::confint.glm() tied to eval()

Hi all,

Has anyone had a chance to look at this and either validate my finding
or tell me that my brain has turned to mush?

Either would be welcome...  :-)
Looks like MASS:::profile.glm could be simplified to

profile.glm <- function(fitted, which = 1:p, alpha = 0.01,
                         maxsteps = 10, del = zmax/5, trace = FALSE, ...)
{
     Pnames <- names(B0 <- coefficients(fitted))
     pv0 <- t(as.matrix(B0))
     p <- length(Pnames)
     if(is.character(which)) which <- match(which, Pnames)
     summ <- summary(fitted)
     std.err <- summ$coefficients[, "Std. Error"]
     mf <- model.frame(fitted)

and this will then work.  That used not to work, including when the 
function was last updated.  The reason it works is that model=TRUE is now 
the default on all the model-fitting functions, as re-creating the model 
frame proved to be too error-prone.
Thanks,

Marc

On Wed, 2006-12-13 at 13:53 -0600, Marc Schwartz wrote:
Greetings all,

I was in the process of creating a function to generate profile
likelihood confidence intervals for a proportion using a binomial glm.
This is a component of a larger function to generate and plot confidence
intervals for proportions using the above, along with bootstrap (BCa),
Wilson and Exact to visually demonstrate the variation across the
methods to some folks.

I had initially tried this using:

  R version 2.4.0 Patched (2006-11-19 r39944)

However, I just verified that it still occurs in:

  R version 2.4.1 RC (2006-12-11 r40157)

In both cases, I started R using '--vanilla' and this is under FC6,
fully updated.

Using the following code, it runs fine:

x <- 14
n <- 35
conf <- 0.95
DF <- data.frame(y = x / n, weights = n)
mod <- glm(y ~ 1, weights = weights, family = binomial(), data = DF)

pl.ci <- plogis(confint(mod, level = conf))
Waiting for profiling to be done...

pl.ci
    2.5 %    97.5 %
0.2490412 0.5651094

However, when I encapsulate the above in a function:

PropCI <- function(x, n, conf = 0.95)
{
  DF <- data.frame(y = x / n, weights = n)
  mod <- glm(y ~ 1, weights = weights, family = binomial(), data = DF)
  plogis(confint(mod, level = conf))
}

I get the following:

# Nothing else in the current global environment
ls()
[1] "PropCI"

PropCI(14, 35)
Waiting for profiling to be done...
Error in inherits(x, "data.frame") : object "DF" not found

If however, I create a 'DF' in the global environment (which may NOT be
the same 'DF' values as within the function!!):

DF <- data.frame(y = 14 / 35, weights = 35)

DF
    y weights
1 0.4      35

ls()
[1] "DF"     "PropCI"

PropCI(14, 35)
Waiting for profiling to be done...
    2.5 %    97.5 %
0.2490412 0.5651094

To my point above about the 'DF' in the global environment not being the
same as the 'DF' within the function body:

DF <- data.frame(y = 5 / 40, weights = 40)
DF
      y weights
1 0.125      40

PropCI(14, 35)
Waiting for profiling to be done...
    2.5 %    97.5 %
0.3752233 0.4217556

Rather scary I think....

So, this suggested a problem in where confint.glm() was looking for
'DF'.

I subsequently traced through the code using debug(), having removed
'DF' from the global environment:

debug(MASS:::confint.glm)
PropCI(14, 35)
debugging in: confint.glm(mod, level = conf)

...

debug: object <- profile(object, which = parm, alpha = (1 - level)/4,
    trace = trace)
Browse[1]>
Error in inherits(x, "data.frame") : object "DF" not found

which brought me to profile.glm():

debug(MASS:::profile.glm)
PropCI(14, 35)
Waiting for profiling to be done...
debugging in: profile.glm(object, which = parm, alpha = (1 - level)/4,
trace = trace)

...

debug: mf <- update(fitted, method = "model.frame")
Browse[1]>
Error in inherits(x, "data.frame") : object "DF" not found

which brought me to update.default():

debug(update.default)
PropCI(14, 35)
Waiting for profiling to be done...
debugging in: update.default(fitted, method = "model.frame")

...

debug: if (evaluate) eval(call, parent.frame()) else call
Browse[1]>
Error in inherits(x, "data.frame") : object "DF" not found

which brought me to eval():

debug(eval)
PropCI(14, 35)
debugging in: eval(mf, parent.frame())

...

debugging in: eval(mf, parent.frame())
debug: .Internal(eval(expr, envir, enclos))
Browse[1]>
Error in inherits(x, "data.frame") : object "DF" not found

Unfortunately, at this point, largely due to lack of sleep of late, my
eyes are sufficiently glazed over and my brain sufficiently vapor locked
that the resolution is not immediately clear and I have not yet dug into
the C code for eval().

Presumably, either I am missing something fundamental here, or there is
a bug where, environment-wise, these respective functions are looking in
the wrong place for 'DF', probably based upon the default environment
arguments noted above.

Any comments from a fresh pair of eyes would be most welcome.

Regards and thanks,

Marc Schwartz

______________________________________________
R-devel at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595

Curious finding in MASS:::confint.glm() tied to eval()

Thread (5 messages)