[R] Getting AIC from lrm in Design package

[Kyle Werner]

Dear David,

Thank you for the reference to Frank Harrell's excellent text. I will
read up to correct my statistical deficiencies offline.

Thank you.


On Sun, Oct 25, 2009 at 1:24 PM, David Winsemius <dwinsemius at comcast.net> wrote:
>
> 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
>
>