[R] Firths bias correction for log-linear models

[sdfrost]

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

__
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.html

Re: [R] Firths bias correction for log-linear models

[David Firth]

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

__
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.html