R help - Jan 2006 - Firths bias correction for log-linear models

Dear R-Help List,

I'm trying to implement Firth's (1993) bias correction for log-linear
models.
Firth (1993) states that such a correction can be implemented by supplementing
the data with a function of h_i, the diagonals from the hat matrix, but
doesn't
provide further details. I can see that for a saturated log-linear model,  h_i=1
for all i, hence one just adds 1/2 to each count, which is equivalent to the
Jeffrey's prior, but I'd also like to get bias corrected estimates for
other log
linear models. It appears that I need to iterate using GLM, with the weights
option and h_i, which I can get from the function hatvalues. For logistic
regression, this can be performed by splitting up each observation into response
and nonresponse, and using weights as described in Heinze, G. and Schemper, M.
(2002), but I'm unsure of how to implement the analogue for log-linear
models. A
procedure using IWLS is described by Firth (1992) in Dodge and Whittaker (1992),
but this book isn't in the local library, and its $141+ on Amazon. I've
tried
looking at the code in the logistf and brlr libraries, but I haven't had any
(successful) ideas. Can anyone help me in describing how to implement this in R?

Thanks!
Simon

On 12 Jan 2006, at 20:54, sdfrost at ucsd.edu wrote:
> Dear R-Help List,
>
> I'm trying to implement Firth's (1993) bias correction for
log-linear
> models.
> Firth (1993) states that such a correction can be implemented by 
> supplementing
> the data with a function of h_i, the diagonals from the hat matrix, 
> but doesn't
> provide further details. I can see that for a saturated log-linear 
> model,  h_i=1
> for all i, hence one just adds 1/2 to each count, which is equivalent 
> to the
> Jeffrey's prior, but I'd also like to get bias corrected estimates
for
> other log
> linear models. It appears that I need to iterate using GLM, with the 
> weights
> option and h_i, which I can get from the function hatvalues. For 
> logistic
> regression, this can be performed by splitting up each observation 
> into response
> and nonresponse, and using weights as described in Heinze, G. and 
> Schemper, M.
> (2002), but I'm unsure of how to implement the analogue for log-linear 
> models. A
> procedure using IWLS is described by Firth (1992) in Dodge and 
> Whittaker (1992),
> but this book isn't in the local library, and its $141+ on Amazon. 
> I've tried
> looking at the code in the logistf and brlr libraries, but I haven't 
> had any
> (successful) ideas. Can anyone help me in describing how to implement 
> this in R?
I don't recommend the adjusted IWLS approach in practice, because that 
algorithm is only first-order convergent.  It is mainly of theoretical 
interest.

The brlr function (in the brlr package) provides a template for a more 
direct approach in practice.  The central operation there is an 
application of optim(), with objective function
  - (l + (0.5 * log(detinfo)))
in which l is the log likelihood and detinfo is the determinant of the 
Fisher information matrix.  In the case of a Poisson log-linear model, 
the Fisher information is, using standard GLM-type notation, t(X) %*% 
diag(mu) %*% X.  It is straightforward to differentiate this penalized 
log-likelihood function, so (as in brlr) derivatives can be supplied 
for use with a second-order convergent optim() algorithm such as BFGS 
(see ?optim for a reference on the algorithm).

I hope that helps.  Please feel free to contact me off the list if 
anything is unclear.

Kind regards,
David

--
Professor David Firth
http://www.warwick.ac.uk/go/dfirth