2

I could use some help wrapping a loop around a function.

I have 26 different models for different plant species that all have identical explanatory variables and structure. Ultimately I want to extract the model coefficients to a table.

First, I created a function to extract the coefficients from one model and put them in a row of an empty dataframe called mod.out. I can run this function for a single model by typing in the model name and a unique row number. 

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)
}

What I would like to do now is write a loop to go through this function for every model to put each set of coefficients in a new row in the mod.out dataframe. The models are glmers. I created a list of all the model names:

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)

I thought I would be able to write a loop function pretty easily to go through it, but I cannot get it to work. I've tried lots of different flavors of the get() and paste() command, but I'm doing something wrong. I think that the problem is in how I am specifying the model name when the function is inside the loops, but I can't figure it out. Any help would be greatly appreciated. Right now I have:

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

I know that there are packages that do things similar to this, but I am working hard to learn functions and loops, so I would really like to do it this way if possible. Thanks!

Skye G.
  • 21
  • 2
  • 2
    Put the models in a list, `mod.name <- list(abam.mort, etc)` then `lapply(seq_along(mod.name), function(i) coefs(mod.name[[i]], i))`. – Rui Barradas Jul 02 '19 at 20:59
  • 1
    It's easier to help you if you include a simple [reproducible example](https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example) with sample input and desired output that can be used to test and verify possible solutions. Note that `c()` doesn't always create a list the way you want it do. It's often better to use `list()` explicitly. And you can't `deparse(substitute())` through multiple variable levels. Better to explicitly track the model names separately. – MrFlick Jul 02 '19 at 21:11
  • @RuiBarradas That code does run without giving me an error (thanks!) but it does not fill the mod.out table out as I was hoping to do. It produces 26 different tables in which only one line is filled with values and the rest are all NAs. – Skye G. Jul 02 '19 at 21:13

1 Answers1

1

As an alternative you might consider this kind of approach:

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

(And a bunch of factor/character conversion warnings.)

Note that this approach relies on a named list of glmer models.

If you stick with your current approach I would at least strongly recommend against extracting model coefficients directly from model objects using $; the package authors can change the internal structure of those objects rendering your code broken. You can get the (fixed) coefficients in this case with something like coef(summary(abam.mort)).

joran
  • 169,992
  • 32
  • 429
  • 468