CODEWITHSUNDEEP
rregressionextractordinalmass-package
theforestecologist
theforestecologist

Reputation: 4937

How to extract AIC from polr summary output in R

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.

Upvotes: 1

Views: 1106

Answers (2)

rr19
rr19

Reputation: 111

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.

Upvotes: -1

Rui Barradas
Rui Barradas

Reputation: 76460

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

Upvotes: 3

Related Questions