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

From: Marc Schwartz <marc_schwartz_at_comcast.net>
Date: Mon 18 Dec 2006 - 14:10:50 GMT

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



R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel Received on Tue Dec 19 08:37:39 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 Mon 18 Dec 2006 - 22:30:56 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.