[R] Getting AIC from lrm in Design package

[David Winsemius]

On Oct 25, 2009, at 12:55 PM, Kyle Werner wrote:

> David,
>
> Thank you for your reply. I am not using glm, but instead lrm.

Does not matter. "lrm" is giving you the same output as would glm with  
a logistic link.

> I am
> consulting the documentation to try to parse out what the output
> "Model L.R." actually means:
> http://lib.stat.cmu.edu/S/Harrell/help/Design/html/lrm.fit.html
> ("model likelihood ratio chi-square")
>
> From my read of the documentation, it appears that the "Model L.R."
> output by lrm is not the deviance, but
>
> 2 * ln( [model deviance] / [null deviance] )

Believe it to be 2 * ln( [null deviance] / [model deviance] ), which  
is why you and I differ as to sign. On page 183 Harrell formulates the  
L.R stat as:

-2log(L at H0/ Lat MLEs)

>
> This is why I believe I am correct to be talking about likelihood
> ratios instead of deviance, but I am unsure of this - which really
> relates back to the core of my question. Whether or not I have the
> sign wrong in the AIC formulation is dependent upon whether or not the
> - is incorporated into the calculation of the "Model L.R.", above,
> which is one part of my original question.
>
> The AIC formula is, generally, AIC = -2*ln(likelihood ratio) + 2k,
> with the best model (assuming the same observations) having the lowest
> AIC. I hope that my understanding of this fundamental formula is
> correct, but please let me know if not.

I have offered my opinion as to your misunderstanding, and have  
suggested specific references to what Harrell has written about the  
subject.

>
> Thanks.
>
> On Sun, Oct 25, 2009 at 10:51 AM, David Winsemius
> <dwinsemius at comcast.net> wrote:
>>
>> On Oct 25, 2009, at 9:24 AM, Kyle Werner wrote:
>>
>>> I am trying to obtain the AICc after performing logistic regression
>>> using the Design package. For simplicity, I'll talk about the AIC. I
>>> tried building a model with lrm, and then calculating the AIC as
>>> follows:
>>>
>>> likelihood.ratio <-
>>> unname(lrm(succeeded~var1+var2,data=scenario,x=T,y=T)$stats["Model
>>> L.R."]) #Model likelihood ratio???
>>> model.params <- 2 #Num params in my model
>>> AIC = -2*log(likelihood.ratio) + 2 * model.params
>>
>> You might want to check your terminology. A single model has a  
>> deviance. You
>> construct a likelihood ratio as twice the logged ratio between two
>> likelihoods (deviances) (which then turns into a difference on the  
>> log
>> scale). And don't you have the sign wrong on that expression for  
>> AIC? I
>> thought you penalized (i.e. subtracted) for added degrees of  
>> freedom? (There
>> is an implicit base model, so if you define AIC as a difference  
>> between L1 +
>> 2p1 and L2+2p2 you would get (L1-L2) + (0 -2p2) = (LR - 2p2). See p  
>> 202 of
>> Harrell's "Regression Modeling Strategies".)
>>
>>>
>>> However, this is almost certainly the wrong interpretation. When I
>>> replace var1 and var2 by runif(N,0,1) (that is, random variables), I
>>> obtain better (lower) AICs than when I use real var1 and var2 that  
>>> are
>>> known to be connected with the outcome variable. Indeed, when I use
>>> GLM instead of LRM, the real model has a lower AIC than that with  
>>> the
>>> predictors replaced by random values. Therefore, it appears that  
>>> lrm's
>>> "Model L.R." is not actually the model likelihood ratio, but instead
>>> something else.
>>>
>>> Going to the Design documentation, lrm.fit states that "Model  
>>> L.R." is
>>> the model likelihood-ratio chi-square. Does this mean that it is
>>> returning 2*log(likelihood)? If so, AIC becomes
>>> AIC = -[Model L.R.] + 2*model.params
>>>
>>> Can anyone confirm that this final formula for obtaining the AIC  
>>> from
>>> lrm is correct?
>>
>>
>> I would not confirm it. What sources are you consulting?
>>
>> --
>> David Winsemius, MD
>> Heritage Laboratories
>> West Hartford, CT
>>
>>

David Winsemius, MD
Heritage Laboratories
West Hartford, CT