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

From: Peter Dalgaard <p.dalgaard_at_biostat.ku.dk>
Date: Sat 16 Dec 2006 - 19:35:39 GMT

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... :-)
>
>

Scoping issues can turn anyones brain to mush! This sort of thing easily happens with code ported from S-PLUS which has different scoping rules.

Inside MASS:::confint.glm we have a call to profile which has a call to update. I think one or both of these need to be evaluated in the parent environment.

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



R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel Received on Sun Dec 17 06:37:53 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 Sat 16 Dec 2006 - 20:31:07 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.