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