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

From: Prof Brian Ripley <ripley_at_stats.ox.ac.uk>
Date: Mon 18 Dec 2006 - 12:53:57 GMT

On Sat, 16 Dec 2006, Marc Schwartz wrote:

> 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@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>

-- 
Brian D. Ripley,                  ripley@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

______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Received on Tue Dec 19 00:02:28 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 Mon 18 Dec 2006 - 22:30:56 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.