Re: [R] Pointwise Confidence Bounds on Logistic Regression

From: Rolf Turner <r.turner_at_auckland.ac.nz>
Date: Thu, 19 Jun 2008 10:22:05 +1200

On 19/06/2008, at 9:32 AM, Prof Brian Ripley wrote:

> On Thu, 19 Jun 2008, Rolf Turner wrote:
>
>>
>> On 19/06/2008, at 8:08 AM, Bryan Hanson wrote:
>>
>>> Hi all. I hope I have my terminology right here...
>>> For a simple lm, one can add “pointwise confidence bounds” to a
>>> fitted line
>>> using something like
>>>> predict(results.lm, newdata = something, interval = "confidence")
>>> (I'm following DAAG page 154-155 for this)
>>> I would like to do the same thing for a glm of the logistic
>>> regression type,
>>> for instance, the example in MASS pg 190-192 (available in the
>>> help page for
>>> predict.glm).
>>> However, predict.glm does not have the same kind of features as
>>> "plain old"
>>> predict, i.e. One cannot specify interval = "confidence"
>>
>> I guess that one reason for that is that prediction intervals
>> rarely if ever make sense with generalized linear models. So only
>> one kind of interval is in effect possible.
>>>> From what I've read, "pointwise confidence bounds" are computed
>>>> from the
>>> SE's for each point. However, I don't see quite where to extract
>>> this
>>> information with a glm
>>> So, is there an existing function that does what I am describing
>>> for a glm,
>>> or can someone point me in the right direction to start writing
>>> my own?
>>
>> Use predict(<whatever>,type="response",se.fit=TRUE). You get a
>> list with
>> three components, the first two of which are the fitted values and
>> their
>> standard errors. (The third is the ``scale'' factor, usually/
>> often equal to 1.)
>
> I would suggest rather computing confidence intervals on linear
> predictor scale
> and transforming those to response scale. That way you will not
> get e.g. negative
> values for probabilities in a logistic regression.

I think that for once Brian, you are being too kind! :-) My advice was even dumber than
you indicate.

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.

Doing what I suggested gives C and D (say) such that P(C <= E(expit (eta-hat)) <= D) = 0.95
(where ``E'' means expected value).

But of course, although E(eta-hat) = eta, it is *NOT* true that E (expit(eta-hat)) = expit(eta).
So what I proposed does NOT give the confidence interval that is desired.

Duhhhhh. Apologies (to Bryan Hanson) for being misleading.

        cheers,

                Rolf

[It probably doesn't make a *lot* of practical difference ``usually'' --- and the standard
errors are approximations anyway ... but one might as well get it right. Sigh.]

######################################################################
Attention:
This e-mail message is privileged and confidential. If you are not the intended recipient please delete the message and notify the sender. Any views or opinions presented are solely those of the author.

This e-mail has been scanned and cleared by MailMarshal www.marshalsoftware.com

######################################################################

______________________________________________
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 Wed 18 Jun 2008 - 22:26:24 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 - 15:30:46 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