I would like to figure out the equations for calculating "AIC" in both "step() function" and "AIC () function". They are different. Then I just type "step" in the R console, and found the "AIC" used in "step() function" is "extractAIC". I went to the R help, and found: "The criterion used is AIC = - 2*log L + k * edf, where L is the likelihood and edf the equivalent degrees of freedom (i.e., the number of free parameters for usual parametric models) of fit. For linear models with unknown scale (i.e., for lm and aov), -2log L is computed from the deviance and uses a different additive constant to logLik and hence AIC. If RSS denotes the (weighted) residual sum of squares then extractAIC uses for - 2log L the formulae RSS/s - n (corresponding to Mallows' Cp) in the case of known scale s and n log (RSS/n) for unknown scale. AIC only handles unknown scale and uses the formula n log (RSS/n) - n + n log 2? - sum log w where w are the weights." Now, my question is what code I should use to look at the exact calculation process in the AIC()function and extractAIC() function in R? Thanks! Dana -- View this message in context: http://www.nabble.com/AIC-function-and-Step-function-tp20728043p20728043.html Sent from the R help mailing list archive at Nabble.com.

Hi Dana, Of course the only true way to know what the AIC calculations are is to read the source code. From within R, what with namespaces, that is becoming increasingly difficult. The AIC() function is not too hard to find from R.> AICfunction (object, ..., k = 2) UseMethod("AIC") <environment: namespace:stats> So now we know that AIC is in the stats package.> objects("package:stats", all=TRUE)[1] ".MFclass" ".checkMFClasses" ".getXlevels" "AIC" [5] "ARMAacf" "ARMAtoMA" "Box.test" "C" so all we are seeing is the AIC generic method shown above. Very often an S3 method will have a 'default' method, so let's see if we get lucky:> stats:::AIC.defaultfunction (object, ..., k = 2) { ll <- if ("stats4" %in% loadedNamespaces()) stats4:::logLik else logLik if (length(list(...))) { object <- list(object, ...) val <- lapply(object, ll) val <- as.data.frame(t(sapply(val, function(el) c(attr(el, "df"), AIC(el, k = k))))) names(val) <- c("df", "AIC") Call <- match.call() Call$k <- NULL row.names(val) <- as.character(Call[-1]) val } else AIC(ll(object), k = k) } <environment: namespace:stats> so there is the basic AIC method. There are some references to logLik, so let's look for a 'logLik' method too.> stats:::AIC.logLikfunction (object, ..., k = 2) -2 * c(object) + k * attr(object, "df") <environment: namespace:stats> If anyone knows how to reveal all these methods without having to guess, I'd like to know of easier ways to find these methods. objects() doesn't reveal them, showMethods() doesn't appear to, and in my attempts, this occurred:> methods(class = "AIC")Error : cannot allocate vector of size 1.9 Gb R(352,0xa0350074) malloc: *** mmap(size=2023526400) failed (error code=12) *** error: can't allocate region *** set a breakpoint in malloc_error_break to debug Error in as.environment(pos) : no item called "get(".__S3MethodsTable__.", envir = asNamespace(i))" on the search list R(352,0xa0350074) malloc: *** mmap(size=2023526400) failed (error code=12) *** error: can't allocate region *** set a breakpoint in malloc_error_break to debug>So it's just not easy finding these functions at the command prompt. Thankfully R is open source, so I downloaded the latest R-2.8.0.tar.gz to directory /Volumes/KilroyHD/kilroy/Software/R/R-2.8.0/ and untarred it all. Navigate to /Volumes/KilroyHD/kilroy/Software/R/R-2.8.0/R-2.8.0/src/library/stats/R and you're looking at all the stats package R scripts. There's 155 files in there, one named AIC.R so that looks promising. Look in AIC.R and you'll see the two functions shown above, but not the extractAIC bunch. Where to look for the rest? There are many ways to look for extractAIC - you could grep the files from a unix shell, etc. Here's a script I use to allow easier perusal of R source code. This script will concatenate all the R source files in a directory into one file, with file name and path descriptors between files in the output. ### 8< ------------------------ ### inputDir <- "/Volumes/KilroyHD/kilroy/Software/R/R-2.8.0/R-2.8.0/src/library/stats/R" outputDir <- paste(inputDir, "/temp", sep = "") outputFile <- "AllRScriptFiles.R" outputDirFile <- paste(outputDir, outputFile, sep = "/") if ( !file.exists(outputDir) ) { dir.create(outputDir) } cat(paste("### All files from dir", inputDir, "\n\n"), file = outputDirFile) ## Get list of all files in dir fileNames <- list.files(path = inputDir, pattern = "\\.R$", full.names = TRUE) ofcon <- file(outputDirFile, open = "at") for (fi in seq(along = fileNames)) { fni <- fileNames[fi] fc <- readLines(fni, n = -1) cat(paste("\n\n##### File", fni, " #####\n\n"), file = outputDirFile, append = TRUE) writeLines(fc, con = ofcon) flush(ofcon) } close(ofcon) ### 8< ------------------------ ### After running this script, navigate to /Volumes/KilroyHD/kilroy/Software/R/R-2.8.0/R-2.8.0/src/library/stats/R/temp with your favourite editor, and edit the file AllRScriptFiles.R (Of course you will substitute your own directory names in all this.) Find extractAIC, and you'll see it came from file add.R You'll see extractAIC methods for classes coxph, survreg, glm, ... and you can see all the calculations. extractAIC <- function(fit, scale, k = 2, ...) UseMethod("extractAIC") extractAIC.coxph <- function(fit, scale, k = 2, ...) { ## seems that coxph sometimes gives one and sometimes gives two values ## for loglik edf <- length(fit$coef) loglik <- fit$loglik[length(fit$loglik)] c(edf, -2 * loglik + k * edf) } extractAIC.survreg <- function(fit, scale, k = 2, ...) { edf <- sum(fit$df) c(edf, -2 * fit$loglik[2] + k * edf) } extractAIC.glm <- function(fit, scale = 0, k = 2, ...) { n <- length(fit$residuals) edf <- n - fit$df.residual aic <- fit$aic c(edf, aic + (k-2) * edf) } extractAIC.lm <- function(fit, scale = 0, k = 2, ...) { n <- length(fit$residuals) edf <- n - fit$df.residual RSS <- deviance.lm(fit) dev <- if(scale > 0) RSS/scale - n else n * log(RSS/n) c(edf, dev + k * edf) } extractAIC.aov <- extractAIC.lm extractAIC.negbin <- function(fit, scale, k = 2, ...) { n <- length(fit$residuals) edf <- n - fit$df.residual c(edf, -fit$twologlik + k * edf) } HTH Steven McKinney Statistician Molecular Oncology and Breast Cancer Program British Columbia Cancer Research Centre email: smckinney +at+ bccrc +dot+ ca tel: 604-675-8000 x7561 BCCRC Molecular Oncology 675 West 10th Ave, Floor 4 Vancouver B.C. V5Z 1L3 Canada -----Original Message----- From: r-help-bounces at r-project.org on behalf of Dana77 Sent: Thu 11/27/2008 6:11 PM To: r-help at r-project.org Subject: [R] AIC function and Step function I would like to figure out the equations for calculating "AIC" in both "step() function" and "AIC () function". They are different. Then I just type "step" in the R console, and found the "AIC" used in "step() function" is "extractAIC". I went to the R help, and found: "The criterion used is AIC = - 2*log L + k * edf, where L is the likelihood and edf the equivalent degrees of freedom (i.e., the number of free parameters for usual parametric models) of fit. For linear models with unknown scale (i.e., for lm and aov), -2log L is computed from the deviance and uses a different additive constant to logLik and hence AIC. If RSS denotes the (weighted) residual sum of squares then extractAIC uses for - 2log L the formulae RSS/s - n (corresponding to Mallows' Cp) in the case of known scale s and n log (RSS/n) for unknown scale. AIC only handles unknown scale and uses the formula n log (RSS/n) - n + n log 2? - sum log w where w are the weights." Now, my question is what code I should use to look at the exact calculation process in the AIC()function and extractAIC() function in R? Thanks! Dana -- View this message in context: http://www.nabble.com/AIC-function-and-Step-function-tp20728043p20728043.html Sent from the R help mailing list archive at Nabble.com. ______________________________________________ R-help at 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.

Hi Dana, Many thanks to Christos Hatzis who sent me an offline response, pointing out the new functions that make this much easier than my last suggestions: methods() and getAnywhere()> methods("extractAIC")[1] extractAIC.aov* extractAIC.coxph* extractAIC.glm* extractAIC.lm* extractAIC.negbin* [6] extractAIC.survreg* Non-visible functions are asterisked> getAnywhere("extractAIC.coxph")A single object matching ?extractAIC.coxph? was found It was found in the following places registered S3 method for extractAIC from namespace stats namespace:stats with value function (fit, scale, k = 2, ...) { edf <- length(fit$coef) loglik <- fit$loglik[length(fit$loglik)] c(edf, -2 * loglik + k * edf) } <environment: namespace:stats>>Thank you Christos. That said, one of the advantages of getting the source code is that it has comments that are stripped out when the code is sourced into R e.g. from the command line> getAnywhere(AIC.default)A single object matching ?AIC.default? was found It was found in the following places registered S3 method for AIC from namespace stats namespace:stats with value function (object, ..., k = 2) { ll <- if ("stats4" %in% loadedNamespaces()) stats4:::logLik else logLik if (length(list(...))) { object <- list(object, ...) val <- lapply(object, ll) val <- as.data.frame(t(sapply(val, function(el) c(attr(el, "df"), AIC(el, k = k))))) names(val) <- c("df", "AIC") Call <- match.call() Call$k <- NULL row.names(val) <- as.character(Call[-1]) val } else AIC(ll(object), k = k) } <environment: namespace:stats>>From the source file# File src/library/stats/R/AIC.R # Part of the R package, http://www.R-project.org # # This program is free software; you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation; either version 2 of the License, or # (at your option) any later version. # # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # A copy of the GNU General Public License is available at # http://www.r-project.org/Licenses/ #### Return the object's value of the Akaike Information Criterion #### (or "An Inf.. Crit..") AIC <- function(object, ..., k = 2) UseMethod("AIC") ## AIC for logLik objects AIC.logLik <- function(object, ..., k = 2) -2 * c(object) + k * attr(object, "df") AIC.default <- function(object, ..., k = 2) { ## AIC for various fitted objects --- any for which there's a logLik() method: ll <- if("stats4" %in% loadedNamespaces()) stats4:::logLik else logLik if(length(list(...))) {# several objects: produce data.frame object <- list(object, ...) val <- lapply(object, ll) val <- as.data.frame(t(sapply(val, function(el) c(attr(el, "df"), AIC(el, k = k))))) names(val) <- c("df", "AIC") Call <- match.call() Call$k <- NULL row.names(val) <- as.character(Call[-1]) val } else AIC(ll(object), k = k) } Steven McKinney Statistician Molecular Oncology and Breast Cancer Program British Columbia Cancer Research Centre email: smckinney +at+ bccrc +dot+ ca tel: 604-675-8000 x7561 BCCRC Molecular Oncology 675 West 10th Ave, Floor 4 Vancouver B.C. V5Z 1L3 Canada