##このスクリプトは、大部分は群馬大学の青木先生が作ったスクリプトを参考にさせていただいています。ここにお礼を申し上げます。 ##なので、このスクリプトはご自由にお使いください。ただし、動作補償はいたしかねます。 ##古いRでは動かないことがあるようです。できるだけ最新版のR上でお使いください。 ## ##作成者:飯島勇人(2007/4/14) #glmヴァージョン akaikeW <- function(df,family) { nv <- ncol(df)-1 name <- names(df) n <- 2^nv-1 str <- character(nv) result1 <- numeric(n) result2 <- character(n) for (i in 1:n) { k <- i m <- 0 for (j in 1:nv) { if (k%%2) { m <- m+1 str[m] <- name[j] } k <- k %/% 2 } form <- reformulate(str[1:m], name[nv+1]) result1[i] <- AIC(result <- glm(form, family)) result2[i] <- paste(c(name[nv+1], "~", paste(str[1:m], collapse="+")), collapse="") } form2 <- reformulate(1, name[nv+1]) nullAIC <- AIC(result <- glm(form2, family)) nullFormula <- paste(name[nv+1], "~ 1") result1 <- c(result1, nullAIC) result2 <- c(result2, nullFormula) o <- order(result1) table0 <- cbind("Formula"=result2[o], "AIC"=result1[o]) table1 <- as.data.frame(table0) table1$AIC <- as.character(table1$AIC) table1$AIC <- as.numeric(table1$AIC) table1$DeltaAIC <- table1$AIC - min(table1$AIC) table1$w <- exp(-0.5*table1$DeltaAIC)/sum(exp(-0.5*table1$DeltaAIC)) for (p in 1:nv) { table1[,4+p] <- regexpr(name[p], table1$Formula) table1[,4+p] <- ifelse(table1[,4+p]==-1, 0, 1) table1[,4+p] <- table1[,4+p]*table1$w } names(table1) <- c("Formula", "AIC", "DeltaAIC", "w", name[1:nv]) RIV <- as.numeric(nv) for (q in 1:nv) { RIV[q] <- sum(table1[,4+q]) } RIV <- as.data.frame(RIV) rlabel <- c(name[1:nv]) rownames(RIV) <- rlabel final <- list(table1, RIV) name <- c("Value", "RIV") names(final) <- name return(invisible(final)) } ##使うときは #df <- read.csv("data.csv") #dfはデータフレームで、一番右の列に従属変数をおいて下さい。 #attach(df) #上記の関数をそのままコピペ #result <- akaikeW(df,poisson(log)) #familyはpoisson(log)のように。 #result ##のような感じで。 #glm.nbヴァージョン akaikeWnb <- function(df) { nv <- ncol(df)-1 name <- names(df) n <- 2^nv-1 str <- character(nv) result1 <- numeric(n) result2 <- character(n) for (i in 1:n) { k <- i m <- 0 for (j in 1:nv) { if (k%%2) { m <- m+1 str[m] <- name[j] } k <- k %/% 2 } form <- reformulate(str[1:m], name[nv+1]) result1[i] <- AIC(result <- glm.nb(form)) result2[i] <- paste(c(name[nv+1], "~", paste(str[1:m], collapse="+")), collapse="") } form2 <- reformulate(1, name[nv+1]) nullAIC <- AIC(result <- glm.nb(form2)) nullFormula <- paste(name[nv+1], "~ 1") result1 <- c(result1, nullAIC) result2 <- c(result2, nullFormula) o <- order(result1) table0 <- cbind("Formula"=result2[o], "AIC"=result1[o]) table1 <- as.data.frame(table0) table1$AIC <- as.character(table1$AIC) table1$AIC <- as.numeric(table1$AIC) table1$DeltaAIC <- table1$AIC - min(table1$AIC) table1$w <- exp(-0.5*table1$DeltaAIC)/sum(exp(-0.5*table1$DeltaAIC)) for (p in 1:nv) { table1[,4+p] <- regexpr(name[p], table1$Formula) table1[,4+p] <- ifelse(table1[,4+p]==-1, 0, 1) table1[,4+p] <- table1[,4+p]*table1$w } names(table1) <- c("Formula", "AIC", "DeltaAIC", "w", name[1:nv]) RIV <- as.numeric(nv) for (q in 1:nv) { RIV[q] <- sum(table1[,4+q]) } RIV <- as.data.frame(RIV) rlabel <- c(name[1:nv]) rownames(RIV) <- rlabel final <- list(table1, RIV) name <- c("Value", "RIV") names(final) <- name return(invisible(final)) }