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 ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
