[R] Firths bias correction for log-linear models

[David Firth]

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