Firths bias correction for log-linear models

classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

Firths bias correction for log-linear models

Simon Frost
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?

Thanks!
Simon

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Reply | Threaded
Open this post in threaded view
|

Re: Firths bias correction for log-linear models

David Firth-2
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-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html