On 12 Jan 2006, at 20:54,

[hidden email] 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______________________________________________

[hidden email] mailing list

https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide!

http://www.R-project.org/posting-guide.html