From: David Firth <d.firth_at_warwick.ac.uk>

Date: Fri 13 Jan 2006 - 09:00:23 EST

Date: Fri 13 Jan 2006 - 09:00:23 EST

On 12 Jan 2006, at 20:54, sdfrost@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 ______________________________________________ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.htmlReceived on Fri Jan 13 09:07:09 2006

*
This archive was generated by hypermail 2.1.8
: Fri 13 Jan 2006 - 14:19:48 EST
*