Перебирать список моделей и сравнивать соответствие модели с помощью AIC, BIC

У меня есть набор данных с несколькими переменными результата, которые я хотел бы протестировать на основе одного предиктора. При исследовательском анализе я заметил, что некоторые зависимости являются полиномиальными до степени 2, а не линейными. Я хотел бы взглянуть на BIC и AIC, чтобы принять решение о том, какую модель лучше всего использовать.

У меня есть функция lapply, с помощью которой я могу перебирать несколько переменных результата, но теперь я хотел бы добавить вторую модель и сравнить ее соответствие. Однако когда я запускаю эту функцию, она сохраняет только вторую модель, и я не знаю, как добраться до «выходов» для запуска следующей функции. Могу ли я использовать это в одной функции или мне нужны две?

Вот пример из набора данных радужной оболочки глаза

data(iris)
vars <- names(iris[2:4])

models2 <- lapply(vars, function(x) {
  model_list=list(
    mod1=lm(substitute(i ~ Sepal.Length, list(i=as.name(x))), data=iris),
    mod2=lm(substitute(i ~ poly(Sepal.Length,2), list(i=as.name(x))), data=iris))
})

y <- lapply(models2, summary) #This only saves results from mod2

Как мне затем сравнить mod1 с mod2 и извлечь следующую переменную?

data.frame(
  do.call(merge, list(BIC(mod1, mod2), AIC(mod1, mod2))), 
  logLik=sapply(list(mod1, mod2), logLik), 
  anova(mod1, mod2, test='Chisq'))

Если у вас есть модели lm в списке, проверьте compare_performance() в library(performance). Выполняет большую часть тяжелой работы за одну строку.

— 21.05.2024 08:43

r lapply model-fitting

21.05.2024 07:42

Стоит ли изучать PHP в 2026-2027 годах?

Привет всем, сегодня я хочу высказать свои соображения по поводу вопроса, который я уже много раз получал в своем сообществе: "Стоит ли изучать PHP в...

Поведение ключевого слова "this" в стрелочной функции в сравнении с нормальной функцией

В JavaScript одним из самых запутанных понятий является поведение ключевого слова "this" в стрелочной и обычной функциях.

Приемы CSS-макетирования - floats и Flexbox

Здравствуйте, друзья-студенты! Готовы совершенствовать свои навыки веб-дизайна? Сегодня в нашем путешествии мы рассмотрим приемы CSS-верстки - в...

Тестирование функциональных ngrx-эффектов в Angular 16 с помощью Jest

В системе управления состояниями ngrx, совместимой с Angular 16, появились функциональные эффекты. Это здорово и делает код определенно легче для...

Концепция локализации и ее применение в приложениях React ⚡️

Локализация - это процесс адаптации приложения к различным языкам и культурным требованиям. Это позволяет пользователям получить опыт, соответствующий...

Пользовательский скаляр GraphQL

Листовые узлы системы типов GraphQL называются скалярами. Достигнув скалярного типа, невозможно спуститься дальше по иерархии типов. Скалярный тип...

210

Перейти к ответу Данный вопрос помечен как решенный

Ответы 4

Ваш вопрос на самом деле о том, как перебирать списки.

summaries <- sapply(models2, lapply, summary)
colnames(summaries) <- vars
rownames(summaries) <- c("degree = 1", "degree = 2")
#example:
summaries["degree = 2", "Petal.Length"]
#shows the summary

lapply(models2, \(x) {
  res <- AIC(x[[1]], x[[2]])
  rownames(res) <- c("degree = 1", "degree = 2")
  res
  })

#[[1]]
#           df      AIC
#degree = 1  3 179.4644
#degree = 2  4 177.7682
#
#[[2]]
#           df      AIC
#degree = 1  3 387.1350
#degree = 2  4 377.8223
#
#[[3]]
#           df      AIC
#degree = 1  3 183.3711
#degree = 2  4 170.3331

21.05.2024 14:01

Что-то вроде этого?
Вычислите статистику одну за другой и cbind их. Затем перенаправьте их в do.call/rbind, чтобы поместить результаты в один фрейм данных.

data(iris)
vars <- names(iris[2:4])

models2 <- lapply(vars, function(x) {
  model_list=list(
    mod1<-lm(substitute(i ~ Sepal.Length , list(i = as.name(x))), data = iris),
    mod2<-lm(substitute(i ~ poly(Sepal.Length,2), list(i = as.name(x))), data = iris))
})

results <- lapply(models2, \(m) {
  mod1 <- m[[1L]]
  mod2 <- m[[2L]]
  AIC <- AIC(mod1, mod2)
  BIC <- BIC(mod1, mod2)[2L]
  logLik <- sapply(m, logLik)
  anova <- anova(mod1, mod2, test = 'Chisq')
  cbind(AIC, BIC, logLik, anova)
}) |>
  do.call(rbind, args = _)

results
#>       df      AIC      BIC     logLik Res.Df       RSS Df Sum of Sq
#> mod1   3 179.4644 188.4963  -86.73221    148  27.91566 NA        NA
#> mod2   4 177.7682 189.8107  -84.88410    147  27.23618  1 0.6794752
#> mod11  3 387.1350 396.1669 -190.56750    148 111.45916 NA        NA
#> mod21  4 377.8223 389.8649 -184.91116    147 103.36231  1 8.0968484
#> mod12  3 183.3711 192.4030  -88.68553    148  28.65225 NA        NA
#> mod22  4 170.3331 182.3757  -81.16656    147  25.91907  1 2.7331790
#>           Pr(>Chi)
#> mod1            NA
#> mod2  5.549049e-02
#> mod11           NA
#> mod21 6.902974e-04
#> mod12           NA
#> mod22 8.245189e-05

^{Created on 2024-05-21 with reprex v2.1.0}

Что означают эти значения p?

— 21.05.2024 15:00

@UseR10085 UseR10085 Тест anovas на важность члена 2-го порядка при наличии точки пересечения и члена 1-го порядка.

— 21.05.2024 16:45

Было бы здорово, если бы вместо mod были имена переменных...

— 21.05.2024 20:08

@ UseR10085 Почему? AIC, BIC, логарифмическое правдоподобие и ANOVA предназначены для моделей, а не для переменных по отдельности.

— 21.05.2024 22:34

21.05.2024 14:10

Ответ принят как подходящий

Во-первых, сделайте lm красивее, используя do.call и reformulate. Затем lapply над моделями вот так:

> models2 <- lapply(setNames(vars, vars), function(x) {
+   list(
+     mod1=do.call('lm', list(reformulate('Sepal.Length', x), quote(iris))),
+     mod2=do.call('lm', list(reformulate('poly(Sepal.Length, 2)', x), quote(iris)))
+   )
+ })
> 
> (res <- lapply(models2, \(x) data.frame(
+   with(x, do.call('merge', list(BIC(mod1, mod2), AIC(mod1, mod2)))),
+   logLik=with(x, sapply(list(mod1, mod2), logLik)),
+   with(x, anova(mod1, mod2))
+ )))
$Sepal.Width
  df      BIC      AIC    logLik Res.Df      RSS Df Sum.of.Sq        F     Pr..F.
1  3 188.4963 179.4644 -86.73221    148 27.91566 NA        NA       NA         NA
2  4 189.8107 177.7682 -84.88410    147 27.23618  1 0.6794752 3.667285 0.05743267

$Petal.Length
  df      BIC      AIC    logLik Res.Df      RSS Df Sum.of.Sq        F      Pr..F.
1  3 396.1669 387.1350 -190.5675    148 111.4592 NA        NA       NA          NA
2  4 389.8649 377.8223 -184.9112    147 103.3623  1  8.096848 11.51519 0.000887343

$Petal.Width
  df      BIC      AIC    logLik Res.Df      RSS Df Sum.of.Sq        F      Pr..F.
1  3 192.4030 183.3711 -88.68553    148 28.65225 NA        NA       NA          NA
2  4 182.3757 170.3331 -81.16656    147 25.91907  1  2.733179 15.50122 0.000126881

Можно дополнительно rbind.

> do.call('rbind', res)
               df      BIC      AIC     logLik Res.Df       RSS Df Sum.of.Sq         F      Pr..F.
Sepal.Width.1   3 188.4963 179.4644  -86.73221    148  27.91566 NA        NA        NA          NA
Sepal.Width.2   4 189.8107 177.7682  -84.88410    147  27.23618  1 0.6794752  3.667285 0.057432671
Petal.Length.1  3 396.1669 387.1350 -190.56750    148 111.45916 NA        NA        NA          NA
Petal.Length.2  4 389.8649 377.8223 -184.91116    147 103.36231  1 8.0968484 11.515191 0.000887343
Petal.Width.1   3 192.4030 183.3711  -88.68553    148  28.65225 NA        NA        NA          NA
Petal.Width.2   4 182.3757 170.3331  -81.16656    147  25.91907  1 2.7331790 15.501223 0.000126881

Данные:

> data(iris)
> vars <- names(iris[2:4])

21.05.2024 16:59

В вопросе models2 указан список списков, поэтому summary пытается получить сводку списков, а не сводку объектов lm. Используя models2, мы можем переписать код, который его использует, вот так. Обратите внимание, что broom::glance создает однострочный тиббл с многочисленными статистическими данными, и, изменяя stats, мы можем выбрать, какие из них мы хотим сохранить.

library(broom)
library(dplyr)
library(purrr)

stats <- c("AIC", "BIC", "logLik")
pval <- function(x, y) anova(x, y, test = "Chisq")$"Pr(>Chi)"[2]

models2 %>%
 setNames(vars) %>%
  map_dfr(with, c(m1 = glance(mod1)[stats], m2 = glance(mod2)[stats],
    pval = pval(mod1, mod2)), .id = "response")

предоставление

# A tibble: 3 × 8
  response     m1.AIC m1.BIC m1.logLik m2.AIC m2.BIC m2.logLik      pval
  <chr>         <dbl>  <dbl>     <dbl>  <dbl>  <dbl>     <dbl>     <dbl>
1 Sepal.Width    179.   188.     -86.7   178.   190.     -84.9 0.0555   
2 Petal.Length   387.   396.    -191.    378.   390.    -185.  0.000690 
3 Petal.Width    183.   192.     -88.7   170.   182.     -81.2 0.0000825

Вот такую статистику производит glance. Здесь статистика означает F-статистику.

names(glance(lm(demand ~., BOD)))
##  [1] "r.squared"     "adj.r.squared" "sigma"         "statistic"    
##  [5] "p.value"       "df"            "logLik"        "AIC"          
##  [9] "BIC"           "deviance"      "df.residual"   "nobs"

Обратите внимание, что модели2 в этом вопросе хороши, но можно было бы написать так. Сначала мы определяем fm, что такое mod1, а затем просто обновляем его, чтобы получить mod2.

models2 <- lapply(vars, function(x) {
  fm <- do.call("lm", list(reformulate("Sepal.Length", x), quote(iris)))
  list(mod1 = fm,
       mod2 = update(fm, . ~ poly(Sepal.Length, 2)))
})

Спасибо за предложение кода models2, это намного более чистый код. Ваш код для получения сводки AIC и т. д. тоже работал хорошо.

— 23.05.2024 03:48

21.05.2024 20:53