On 12 Jan 2006, at 20:54, [EMAIL PROTECTED] 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 ______________________________________________ [email protected] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
