Обертывание цикла вокруг функции - PullRequest
Новая
       5

Обертывание цикла вокруг функции

2 голосов
/ 02 июля 2019

Я мог бы использовать некоторую помощь, оборачивая цикл вокруг функции.

У меня есть 26 различных моделей для разных видов растений, которые имеют одинаковые объясняющие переменные и структуру. В конечном итоге я хочу извлечь коэффициенты модели в таблицу.

Сначала я создал функцию для извлечения коэффициентов из одной модели и поместил их в ряд пустого фрейма данных с именем mod.out. Я могу запустить эту функцию для отдельной модели, введя название модели и уникальный номер строки. 

coefs<- function(model, row.num){
mod.out[row.num,1]<-strtrim(deparse(substitute(model)), 4)
mod.out[row.num, 2:4]<-summary(model)$coefficients[1, c(1,2,4)] 
mod.out[row.num,5:7]<-summary(model)$coefficients[2, c(1,2,4)]
mod.out[row.num,8:10]<-summary(model)$coefficients[3, c(1,2,4)]
mod.out[row.num,11:13]<-summary(model)$coefficients[4, c(1,2,4)]
mod.out[row.num,14]<-summary(model)$optinfo$val[1]
return(mod.out)
}

Что я хотел бы сделать сейчас, так это написать цикл, чтобы пройти через эту функцию для каждой модели, чтобы поместить каждый набор коэффициентов в новую строку в кадре данных mod.out. Модели блестят. Я создал список всех названий моделей: 

mod.name<-c(abam.mort, abco.mort, abgr.mort, abla.mort, acma.mort, arme.mort, cade.mort, chch.mort, chla.mort, juoc.mort, laoc.mort, lide.mort, pial.mort, piat.mort, pico.mort, pien.mort, pije.mort, pila.mort, pimo.mort, pipo.mort, psme.mort, quch.mort, thpl.mort, tshe.mort, tsme.mort, umca.mort)

Я думал, что смогу довольно легко написать функцию цикла, чтобы пройти через нее, но я не могу заставить ее работать. Я пробовал много разных разновидностей команд get () и paste (), но я делаю что-то не так. Я думаю, что проблема в том, как я определяю имя модели, когда функция находится внутри циклов, но я не могу понять это. Любая помощь будет принята с благодарностью. Прямо сейчас у меня есть: 

for(i in 1:nrow(mod.out)){
  coefs(mod.name[i], i)}

Я знаю, что есть пакеты, которые делают подобные вещи, но я усердно работаю над изучением функций и циклов, поэтому я действительно хотел бы сделать это таким образом, если это возможно. Спасибо!

1 Ответ

1 голос
/ 03 июля 2019

В качестве альтернативы вы можете рассмотреть такой подход: 

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

gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
                         data = cbpp, family = binomial)
l <- list(mod1 = gm1,mod2 = gm1)

> map_dfr(l,tidy,.id = "model")

# A tibble: 10 x 7
   model term                estimate std.error statistic   p.value group
   <chr> <chr>                  <dbl>     <dbl>     <dbl>     <dbl> <chr>
 1 mod1  (Intercept)          -1.398     0.2312    -6.048  1.468e-9 fixed
 2 mod1  period2              -0.9919    0.3032    -3.272  1.068e-3 fixed
 3 mod1  period3              -1.128     0.3228    -3.495  4.745e-4 fixed
 4 mod1  period4              -1.580     0.4220    -3.743  1.818e-4 fixed
 5 mod1  sd_(Intercept).herd   0.6421   NA         NA     NA        herd 
 6 mod2  (Intercept)          -1.398     0.2312    -6.048  1.468e-9 fixed
 7 mod2  period2              -0.9919    0.3032    -3.272  1.068e-3 fixed
 8 mod2  period3              -1.128     0.3228    -3.495  4.745e-4 fixed
 9 mod2  period4              -1.580     0.4220    -3.743  1.818e-4 fixed
10 mod2  sd_(Intercept).herd   0.6421   NA         NA     NA        herd

(и несколько предупреждений о преобразовании коэффициента / символа.)

Обратите внимание, что этот подход основан на именованном спискемоделей glmer.

Если вы придерживаетесь своего текущего подхода, я бы по крайней мере настоятельно рекомендовал не извлекать коэффициенты модели непосредственно из объектов модели, используя $;авторы пакетов могут изменить внутреннюю структуру этих объектов, делая ваш код поврежденным.Вы можете получить (фиксированные) коэффициенты в этом случае с чем-то вроде coef(summary(abam.mort)).

...