Comparing models with dplyr and broom::glance: How

2019-04-17 12:36发布

问题:

I would like to run each variable in a dataset as a univariate glmer model using the lme4 package in R. I would like to prepare the data with the dplyr/tidyr packages, and organize the results from each model with the broom package (i.e. do(glance(glmer...). I would most appreciate help that stuck within that framework. I'm not that great in R, but was able to produce a dataset that throws an error and has the same structure as the data I'm using:

library(lme4)
library(dplyr)
library(tidyr)
library(broom)

Bird<-c(rep(c(0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0),10))
Stop<-c(rep(seq(1,10), 20))
Count<-c(rep(c(rep(c(1,2), each=10)), each=10))
Route<-c(rep(seq(1,10), each=20))
X1<-rnorm(200, 50, 10)
X2<-rnorm(200, 10, 1)
X3<-c(rep(c(0),200))#trouble maker variable

Data<-data.frame(cbind(Bird, Stop, Count, Route, X1, X2, X3))

Data%>%
  gather(Variable, Value, 5:7)%>%
  group_by(Variable)%>%
  do(glance(glmer(Bird~Value+Stop+(1+Stop|Route/Count), data=.,     family=binomial)))

The last variable produces an error so there is no output. What I would like is it to produce NA values in the output if this occurs, or just skip that variable. I've tried using 'try' to blow past the trouble maker variable:

do(try(glance(glmer(Bird~Value+Stop+(1+Stop|Route/Count), data=.,       family=binomial))))

which it does, but still an output is not produced because it can't coerce a 'try-error' to a data.frame. Unfortunately there is no tryharder function. I've tried some if statements which make sense to me but not the computer. I'm sure I'm not doing it right, but if for example I use:

try(glance(glmer(Bird~Value+Stop+(1+Stop|Route/Count), data=., family=binomial)))->mod
if(is.data.frame(mod)){do(mod)}

I get subscript out of bounds errors. Thanks very much for any input you can provide!

回答1:

Use tryCatch before the call to glance:

zz = Data %>%
  gather(Variable, Value, 5:7) %>%
  group_by(Variable) %>%
  do(aa = tryCatch(glmer(Bird~Value+Stop+(1+Stop|Route/Count), data=., 
                  family=binomial), error = function(e) data.frame(NA))) 


zz %>% 
  glance(aa)