[R] SE of difference in fitted probabilities from logistic model.

[Prof Brian Ripley]

On Fri, 16 May 2008, Rolf Turner wrote:

>
> I am fitting a logistic binomial model of the form
>
> 	glm(y ~ a*x,family=binomial)
>
> where a is a factor (with 5 levels) and x is a continuous predictor.
>
> To assess how much ``impact'' x has, I want to compare the fitted 
> success probability when x = its maximum value with the fitted 
> probability when x = its mean value. (The mean and the max are to be 
> taken by level of the factor ``a'', but that's not really an issue.)
>
> I can of course easily calculate p.hat(x.max) - p.hat(x.mean) using 
> predict() (with type="response").  And I can get the standard error for 
> p.hat(x.max) and p.hat(x.mean) by specifiying se.fit=TRUE. No problem 
> there.
>
> But how can I get a handle on the standard error of the difference?

The same way predict() does?  It used a crude (delta-function) 
transformation of SEs on linear predictor scale.

I am not sure that a difference in probabilites is interesting unless it 
is small compared to the values, when linearization will be most likely be 
satisfactory.  (Why not ratios or log odds, for example?)
Linearization can also be used if the SEs on each probability are small 
but the difference is not.  Otherwise you may as well simulate from the 
joint distribution of the coefficients, compute the difference for the 
similations and summarize.  (That might be as quick to do as to work out 
the analytical approximations.)


> In a linear model this would just be SE(beta_1.hat)*(x.max-x.mean) (where
> beta_1.hat is specific to the particular level of `a' being considered).
> If I am not mistaken.  (Please correct me if I am!)
>
> But in the logistic model, everything is entangled in the inverse link
> function (the ``expit'' function as it is called by some), and I can see
> no way of disentangling.
>
> Is there any way of getting at this?  I figure that simulation/Monte 
> Carlo inference/ parametric bootstrapping would provide a workaround, 
> but before I go that route, can anyone point me to a simpler method? 
> There wouldn't be anything built into R or an R package, would there? 
> (I did a fairly basic RSiteSearch() and came up empty handed.)
>
> Thanks for any tips.
>
> 	cheers,
>
> 		Rolf Turner

-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595