[R] penalized maximum likelihood estimation and logistf

Ravi Varadhan rvaradhan at jhmi.edu
Tue Mar 9 18:22:02 CET 2010

Let us look at a piece of the relevant code from "logistf" package:

        iter <- iter + 1
        XW2 <- crossprod(x, diag(pi * (1 - pi))^0.5)
        Fisher <- crossprod(t(XW2))
        covs <- solve(Fisher)
        H <- crossprod(XW2, covs) %*% XW2
        if (firth) 
            U.star <- crossprod(x, y - pi + diag(H) * (0.5 - 
        else U.star <- crossprod(x, y - pi)
        delta <- as.vector(covs %*% U.star)
        mx <- max(abs(delta))/maxstep
        if (mx > 1) 
            delta <- delta/mx
First, note that this is inside an iterative loop.  The Fisher information matrix is first formed, and then it is inverted to get `covs'.  This is where it can get troublesome.  It is very likely in your case that the `Fisher' information matrix is ill-conditioned, and hence `covs' cannot be computed.  Consequently, the steplength `delta' for the Newton method cannot be determined. 



Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology
School of Medicine
Johns Hopkins University

Ph. (410) 502-2619
email: rvaradhan at jhmi.edu

----- Original Message -----
From: xiaoyan yu <xiaoyan.yu at gmail.com>
Date: Tuesday, March 9, 2010 11:26 am
Subject: [R] penalized maximum likelihood estimation and logistf
To: r-help at r-project.org

> Hi, I got two questions and would really appreciate any help from here.
>  First, is the penalized maximum likelihood estimation(Firth Type Estimation)
>  only fit for binary response (0,1 or TRUE, FALSE)? Can it be applied 
> to
>  multinomial logistic regression?
>  If yes, what's the formula for LL and  U(beta_i)? Can someone point 
> me to
>  the right reference?
>  Second, when I used *logistf *on a dataset with binary response, I 
> got the
>  following error msg:
>  "Error in if (mx > 1) delta <- delta/mx :
>    missing value where TRUE/FALSE needed"
>  What does this msg mean? does it mean the dataset includes missing response?
>  I  checked my dataset several times and didn't find anything other 
> than TRUE
>  or FALSE for the response. Please shed light on this.
>  Thanks,
>  Xiaoyan
>  	[[alternative HTML version deleted]]
>  ______________________________________________
>  R-help at r-project.org mailing list
>  PLEASE do read the posting guide 
>  and provide commented, minimal, self-contained, reproducible code.

More information about the R-help mailing list