Marc Schwartz
2006-Dec-13 19:53 UTC
[Rd] Curious finding in MASS:::confint.glm() tied to eval()
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.ci2.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)> DFy 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) > DFy 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
Marc Schwartz
2006-Dec-16 17:51 UTC
[Rd] Curious finding in MASS:::confint.glm() tied to eval()
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