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