Поиск модели (возвращенной из циклов for) с наименьшим значением AIC в R

Я пытаюсь найти модель с самым низким AIC. Модели возвращаются из двух циклов for, которые делают возможными комбинации столбцов. Я не могу сделать модель возврата функции с самым низким AIC. Код ниже демонстрирует, где я застрял:

rm(list = ls())

data <- iris

data <- data[data$Species %in% c("setosa", "virginica"),]

data$Species = ifelse(data$Species == 'virginica', 0, 1)

mod_headers <- names(data[1:ncol(data)-1])

f <- function(mod_headers){
    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, "Species")
      model <- glm(Species ~., data=data[c(mod_tab_new)], family = binomial(link = "logit"))
    }
    }
  best_model <- model[which(AIC(model)[order(AIC(model))][1])]
  print(best_model)
}

f(mod_headers)

Какие-либо предложения? Спасибо!


r for-loop glm model-comparison
person New2coding    schedule 22.07.2017    source источник
comment
Вы должны прочитать r4ds.had.co.nz/many-models.html   -  person F. Privé    schedule 22.07.2017
comment
rdocumentation.org/packages/MuMIn/versions/1.15.6/ темы/драги   -  person Roland    schedule 22.07.2017

Ответы (3)


arrow_upward
0
arrow_downward

glm() использует итеративный алгоритм наименьших квадратов с повторным взвешиванием. Алгоритм достигает максимального количества итераций, прежде чем сходится — изменение этого параметра помогает в вашем случае:

 glm(Species ~., data=data[mod_tab_new], family = binomial(link = "logit"), control = list(maxit = 50))

Была еще одна проблема с использованием which, я заменил ее на if после того, как каждая модель подошла для сравнения с самым низким AIC на данный момент. Однако я думаю, что есть решения получше, чем этот for-loop подход.

f <- function(mod_headers){
  lowest_aic <- Inf     # added
  best_model <- NULL    # added

  for(i in 1:length(mod_headers)){
    tab <- combn(mod_headers,i)
    for(j in 1:ncol(tab)){
      tab_new <- tab[, j]
      mod_tab_new <- c(tab_new, "Species")
      model <- glm(Species ~., data=data[mod_tab_new], family = binomial(link = "logit"), control = list(maxit = 50))
      if(AIC(model) < lowest_aic){ # added
        lowest_aic <- AIC(model)   # added
        best_model <- model        # added
      }
    }
  }
  return(best_model)
}
person shosaco    schedule 22.07.2017
comment
Это на самом деле работает и кажется наиболее эффективным для поиска моей оптимальной модели с кучей столбцов. Спасибо! - person New2coding; 22.07.2017

arrow_upward
1
arrow_downward

Я заменил ваши циклы for на векторизованные альтернативы

library(tidyverse)
library(iterators)
# Column names you want to use in glm model, saved as list
whichcols <- Reduce("c", map(1:length(mod_headers), ~lapply(iter(combn(mod_headers,.x), by="col"),function(y) c(y))))

# glm model results using selected column names, saved as list
models <- map(1:length(whichcols), ~glm(Species ~., data=data[c(whichcols[[.x]], "Species")], family = binomial(link = "logit")))

# selects model with lowest AIC
best <- models[[which.min(sapply(1:length(models),function(x)AIC(models[[x]])))]]

Выход

Call:  glm(formula = Species ~ ., family = binomial(link = "logit"), 
data = data[c(whichcols[[.x]], "Species")])

Coefficients:
 (Intercept)  Petal.Length  
       55.40        -17.17  

Degrees of Freedom: 99 Total (i.e. Null);  98 Residual
Null Deviance:      138.6 
Residual Deviance: 1.208e-09    AIC: 4
person CPak    schedule 22.07.2017

arrow_upward
0
arrow_downward

Используя свой цикл, просто поместите все модели в один список. Затем вычислите AIC всех этих моделей. Наконец верните модель с минимальным AIC.

f <- function(mod_headers) {

  models <- list()
  k <- 1
  for (i in 1:length(mod_headers)) {
    tab <- combn(mod_headers, i)
    for(j in 1:ncol(tab)) {
      mod_tab_new <- c(tab[, j], "Species")
      models[[k]] <- glm(Species ~ ., data = data[mod_tab_new], 
                         family = binomial(link = "logit"))
      k <- k + 1
    }
  }

  models[[which.min(sapply(models, AIC))]]
}
person F. Privé    schedule 22.07.2017