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

[Marc Schwartz]

On Mon, 2006-12-18 at 12:53 +0000, Prof Brian Ripley wrote:
> On Sat, 16 Dec 2006, 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...  :-)
> 
> Looks like MASS:::profile.glm could be simplified to
> 
> profile.glm <- function(fitted, which = 1:p, alpha = 0.01,
>                          maxsteps = 10, del = zmax/5, trace = FALSE, ...)
> {
>      Pnames <- names(B0 <- coefficients(fitted))
>      pv0 <- t(as.matrix(B0))
>      p <- length(Pnames)
>      if(is.character(which)) which <- match(which, Pnames)
>      summ <- summary(fitted)
>      std.err <- summ$coefficients[, "Std. Error"]
>      mf <- model.frame(fitted)
> 
> and this will then work.  That used not to work, including when the 
> function was last updated.  The reason it works is that model=TRUE is now 
> the default on all the model-fitting functions, as re-creating the model 
> frame proved to be too error-prone.

<snip>

Prof. Ripley,

Thanks for taking time to review this and present a solution.

I created a locally modified version of the VR bundle with this change
in profiles.R, installed it and can confirm that the modification does
work:

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))
}


> plogis(confint(glm(14 / 35 ~ 1, family = binomial, weights = 35)))
Waiting for profiling to be done...
    2.5 %    97.5 % 
0.2490412 0.5651094 

> PropCI(14, 35)
Waiting for profiling to be done...
    2.5 %    97.5 % 
0.2490412 0.5651094 


> plogis(confint(glm(5 / 40 ~ 1, family = binomial, weights = 40)))
Waiting for profiling to be done...
     2.5 %     97.5 % 
0.04672812 0.24963731 

> PropCI(5, 40)
Waiting for profiling to be done...
     2.5 %     97.5 % 
0.04672812 0.24963731 


There was a brief flash in my mind over the weekend, wondering why (in
my examples) 'DF' was being looked for when the model frame could be
used instead. I appreciate your clarification on that point.

Thanks again and regards,

Marc