How to extract AIC from polr summary output in R

Question

I'd like to quickly compare AICs that are provided as output when running summary() on individual polr() models created using the MASS package in R. I have no problem compiling this info, but what I can't figure out is where exactly the AIC info is being stored in the polr model objects themselves.

I've tried using str() and attributes() on my model objects, and I've even tried using getAnywhere("polr") to look at the source code itself. Nothing is standing out to me.

Anyone know how to extract AIC output from summary(polr_mod)?

Example for reference:

library(MASS)
dat <- data.frame(v1 = factor(rep(0:2,each=3),ordered = T), v2 = rep(1:3,each=3))
mod_polr <- polr(v1 ~ v2, data = dat, Hess = T, method = "logistic" )
summary(mod_polr)

Call:
polr(formula = v1 ~ v2, data = dat, Hess = T, method = "logistic")

Coefficients:
     v2 
46.7233 

Intercepts:
      0|1       1|2 
 73.62014 117.05781 

Residual Deviance: 1.560782e-08 
AIC: 6.00

^ See, at the bottom of the output is AIC: 6.00. Where is this stored in the object? I'd like to call/extract it.

Thanks, @danh! That works. Though, I'd love to see an explanation for why that works as a formal answer. — theforestecologist, Oct 21 '21 at 14:00
You could take a look at the guts of `AIC()` here: https://github.com/SurajGupta/r-source/blob/master/src/library/stats/R/AIC.R — Skaqqs, Oct 21 '21 at 14:37

score 3 · Answer 1 · answered Oct 21 '21 at 15:02

AIC is a generic function with no method for objects of class "polr" but with a default method. The default method's code can be seen by running

getAnywhere("AIC.default")

and what it does is to call logLik on its first argument and then to compute the AIC with k = 2.
The number of model parameters is given by

attr(lls, "df")

And the computation is

lls <- logLik(mod_polr)
-2*as.numeric(lls) + k*attr(lls, "df")

In its turn, logLik is also generic but with a method for objects of class "polr".
The code returned by

getAnywhere("logLik.polr")

is very simple, a one-liner. Note that the df in its code was edf in the object "polr". Indented it's the following.

logLik.polr <- function(object, ...) {
  structure(
    -0.5 * object$deviance, 
    df = object$edf, 
    nobs = object[["nobs"]], 
    class = "logLik"
  )
}

The explanation for the deviance, edf and nobs are found in help("polr), section Value.

deviance
the residual deviance.

edf
the (effective) number of degrees of freedom used by the model

And like that section says, nobs is used for stepAIC (but not for logLik.polr).

So this can all be pieced together in one function. Write a AIC method for objects of class "polr" simplifying the code above.

AIC.polr <- function(x, k = 2){
  dev <- x$deviance
  nparams <- x$edf
  dev + k*nparams
}

AIC(mod_polr)
#[1] 6

score 0 · Answer 2 · answered Jan 01 '23 at 21:47

0

I know that tihs is an old question ut I just want to share in case others have the same question and want another option.

Using modelsummary() to create your regression tabel then you also get AIC, BIC and RMSE. If you need example let me know.

answered Jan 01 '23 at 21:47

rr19

79
1
9

How to extract AIC from polr summary output in R

2 Answers2