Подмножество значений хи-квадрат из результатов GAM в R

Я запустил цикл for в R, который генерирует модели для биномиальной GAM для 200 различных комбинаций случайных данных (200 различных значений set.seed).

Цикл for и GAM работают нормально, и они сохраняют модели в соответствующем списке, model[[i]], где каждый элемент списка представляет модель для определенной итерации set.seed.

Я могу запустить summary() для отдельного элемента списка (например, model[[5]]) и получить что-то вроде этого:

Approximate significance of smooth terms:
                      edf Ref.df Chi.sq  p-value    
s(Petal.Width)  1.133e+00      9  5.414 0.008701 ** 
s(Sepal.Length) 8.643e-01      9  2.338 0.067509 .  

---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.361   Deviance explained = 32.7%
-REML = 83.084  Scale est. = 1         n = 160

Поскольку у меня есть 200 таких элементов в моем списке model, мне было интересно, есть ли быстрый способ просмотреть и подсчитать, сколько раз (из общего числа 200) это значение Chi.sq равно 0 для лепестка. Переменная ширины. По сути, поскольку у меня есть GAM, настроенные с bs = "cs", количество раз, когда Chi.sq равно 0, представляет, как часто процесс выбора переменной удалял эту переменную из модели.

Вот очищенная версия кода, который я использовал для цикла for для построения модели:


a <- c(1:200)
model <- list()


for (i in a) {
  
#In my version there's some extra code here that subsets iris based on i 
  
  model[[i]] <- gam(Species~ s(Petal.Width, bs="cs") + 
                      s(Sepal.Length, bs="cs"),
                    data = iris,
                    family = binomial,
                    method = "REML")

}

r subset chi-squared
person pfadenhw    schedule 27.01.2021    source источник
comment
не могли бы вы включить код для вашей модели? Я использовал команду tidy() из пакета broom, которая помещает результаты модели в табличку для сбора вывода, после чего вы можете объединить столбцы.   -  person brian avery    schedule 28.01.2021
comment
@brianavery Сделал некоторые правки — это то, что вы искали?   -  person pfadenhw    schedule 28.01.2021
comment
полезно, спасибо, но какой пакет вы используете для gam -- gam или mgcv?   -  person brian avery    schedule 28.01.2021
comment
@brianavery Я использую mgcv   -  person pfadenhw    schedule 28.01.2021

Ответы (2)


arrow_upward
0
arrow_downward

Я обнаружил, что функция tidy() из пакета broom полезна в подобных ситуациях.
вот ваша модель (я сократил количество итераций до 10, чтобы ускорить работу,
также было выдавать предупреждения с помощью mgcv::gam())

a <- c(1:10)
model <- list()

for (i in a) {
  model[[i]] <- gam(Species~ s(Petal.Width, bs="cs") + 
                      s(Sepal.Length, bs="cs"),
                    data = iris,
                    family = binomial,
                    method = "REML")
}

Если вам нужен или нужен список всех моделей, то этот второй цикл извлечет statistic из каждой модели и объединит их в фрейм данных:

results_df <- data.frame(term=c("Petal.Width","Sepal.Length"))
for(m in model){
  results_df <- bind_cols(results_df, model=tidy(m)$statistic)
}
results_df

с выходом:

          term    model...2    model...3    model...4    model...5    model...6
1  Petal.Width 1.304372e-09 1.304372e-09 1.304372e-09 1.304372e-09 1.304372e-09
2 Sepal.Length 3.844738e-16 3.844738e-16 3.844738e-16 3.844738e-16 3.844738e-16
     model...7    model...8    model...9   model...10   model...11
1 1.304372e-09 1.304372e-09 1.304372e-09 1.304372e-09 1.304372e-09
2 3.844738e-16 3.844738e-16 3.844738e-16 3.844738e-16 3.844738e-16

вы можете переименовать столбцы, чтобы сделать их более красивыми, и т. д.
также вы можете комбинировать моделирование и объединение статистики, если это все, что вы хотите сохранить.

person brian avery    schedule 27.01.2021
comment
Упс - я уменьшил количество переменных, когда упростил код для этого вопроса, поэтому не обращайте внимания на предыдущий комментарий (теперь удаленный). - person pfadenhw; 28.01.2021
comment
Однако теперь у меня есть кадр данных results_df, но статистика в нем не соответствует значениям Chi.sq, которые я получаю, если вручную подтягиваю summary(model[[1]]) - person pfadenhw; 28.01.2021

arrow_upward
0
arrow_downward

Используя информацию от @brian avery и этот вопрос, я придумайте следующее, что точно отвечает на мой первоначальный вопрос. Я уверен, что есть более эффективный способ сделать это с помощью каналов, но это работает для меня.

a <- c(1:200)
model <- list()
sum <- list ()
Chisq <- list()


for (i in a) {
  model[[i]] <- gam(Species~ s(Petal.Width, bs="cs") + 
                      s(Sepal.Length, bs="cs"),
                    data = iris,
                    family = binomial,
                    method = "REML")

  sum[[i]] <- summary(model[[i]])
  Chisq[[i]] <- sum[[i]]$chi.sq

}


results_df <- data.frame(term=c("Petal.Width","Sepal.Length"))

for(i in a){
  results_df <- bind_cols(results_df, Chisq[[i]])
}

rowSums(results_df<0.001)

В последней строке подсчитывается, сколько раз переменные оказывались со значением Chi.sq ниже 0,001. Он отвечает на вопрос, сколько раз (из 200) пакет mgcv автоматически исключал эту переменную из модели.

person pfadenhw    schedule 28.01.2021