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

[Peter Dalgaard]

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