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

From: Marc Schwartz <marc_schwartz_at_comcast.net>
Date: Wed 13 Dec 2006 - 19:53:46 GMT


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@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel Received on Fri Dec 15 13:21:40 2006

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Sun 17 Dec 2006 - 02:31:00 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-devel. Please read the posting guide before posting to the list.