[R] Firths bias correction for log-linear models
David Firth
d.firth at warwick.ac.uk
Thu Jan 12 23:00:23 CET 2006
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
More information about the R-help
mailing list