From: Peter Dalgaard <p.dalgaard_at_biostat.ku.dk>

Date: Sat 16 Dec 2006 - 19:35:39 GMT

R-devel@r-project.org mailing list

https://stat.ethz.ch/mailman/listinfo/r-devel Received on Sun Dec 17 06:37:53 2006

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