0

我正在尝试找到最好的 GLModel（目前仅基于 AIC 分数）。我使用的数据可以在这里访问：

https://drive.google.com/open?id=0B5IgiR_svnKcY25TQ29ZMGN3NVE

下面是一段未能返回最佳模型的代码：

rm(list = ls())

mod_data <- read.csv("mod_data.csv", header = T)

mod_headers <- names(mod_data[3:ncol(mod_data)-1])

f <- function(mod_headers){

  null_model <- glm(newcol ~ 1, data=mod_data, family = binomial(link = "logit"), control = list(maxit = 50))
  best_model <- null_model
  best_aic <- AIC(null_model)

  for(i in 1:length(mod_headers)){
    tab <- combn(mod_headers,i)
    for(j in 1:ncol(tab)){
      tab_new <- c(tab[,j])
      mod_tab_new <- c(tab_new, "newcol")
      model <- glm(newcol ~., data=mod_data[c(mod_tab_new)], family = binomial(link = "logit"), control = list(maxit = 50000))
      if(AIC(model) < best_aic){
        best_model <- model
        best_aic <- AIC(model)
      }
    }
    return(best_model)
  }
}

f(mod_headers)

出于某种原因，它总是卡在只有一个自变量的模型（“var5”，AIC 得分为 678.8），尽管我知道至少有一个更好的模型（“var1”，“var5”和“var8”，AIC 得分为 677.6）返回通过这段代码：

rm(list = ls())

mod_data <- read.csv("mod_data.csv", header = T)

mod_headers <- names(mod_data[3:ncol(mod_data)-1])

library(tidyverse)
library(iterators)

whichcols <- Reduce("c", map(1:length(mod_headers), ~lapply(iter(combn(mod_headers,.x), by="col"),function(y) c(y))))
models <- map(1:length(whichcols), ~glm(newcol ~., data=mod_data[c(whichcols[[.x]], "newcol")], family = binomial(link = "logit")))
print(models[[which.min(sapply(1:length(models),function(x)AIC(models[[x]])))]])

我不知道为什么第一段代码未能返回所需的结果。

感谢您的提示！

4

0 回答 0