Re: [R] Pointwise Confidence Bounds on Logistic Regression

From: Bryan Hanson <hanson_at_depauw.edu>
Date: Thu, 19 Jun 2008 10:42:15 -0400


[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")

R-help_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Thu 19 Jun 2008 - 15:01:51 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Thu 19 Jun 2008 - 18:30:47 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive