[R] Pointwise Confidence Bounds on Logistic Regression

[Bryan Hanson]

[I've ommitted some of the conversation so far...]

> E.g. in a logistic model, with (say) eta = beta_0 + beta_1*x one may
> find, on the
> linear predictor scale, A and B (say) such that P(A <= eta <= B) = 0.95.
> 
> Then P(expit(A) <= expit(eta) <= expit(B)) = 0.95, which is exactly
> what is wanted.

I think I follow the above conceptually, but I don't know how to implement
it, though I fooled around (unsuccessfully) with some of the variations on
predict().

I'm trying to learn this in response to a biology colleague who did
something similar in SigmaPlot. I can already tell that SigmaPlot did a lot
of stuff for him in the background.  The responses are 0/1 of a particular
observation by date.  The following code simulates what's going on (note
that I didn't use 0/1 since this leads to a recognized condition/warning of
fitting 1's and 0's. I've requested Brian's Pattern Recognition book so I
know what the problem is and how to solve it).  My colleague is looking at
two populations in which the "LD50" would differ.  I'd like to be able to
put the "pointwise confidence bounds" on each curve so he can tell if the
populations are really different.

By the way, this code does give a (minor?) error from glm (which you will
see).

Can you make a suggestion about how to get those "confidence bounds" on
there?

Also, is a probit link more appropriate here?

Thanks,  Bryan

x <- c(1:40)
y1 <- c(rep(0.1,10), rep(NA, 10), rep(0.9,20))
y2 <- c(rep(0.1,15), rep(NA, 10), rep(0.9,15))
data <- as.data.frame(cbind(x,y1,y2))
plot(x, y1, pch = 1, ylim = c(0,1), col = "red")
points(x, y2, pch = 3, col = "blue")
abline(h = 0.5, col = "gray")
fit1 <- glm(y1~x, family = binomial(link = "logit"), data, na.action =
na.omit)
fit2 <- glm(y2~x, family = binomial(link = "logit"), data, na.action =
na.omit)
lines(fit1$model$x, fit1$fitted.values, col = "red")
lines(fit2$model$x, fit2$fitted.values, col = "blue")