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