Reputation: 517
I'm running a list of ordered logit models with different variables etc. I want to convert the output into a tidy tibble to use in ggplot etc. (I also want to save the 'regular model output' so I want to do this separately.)
I want to do this in an automated way, using purrr or lapply or some such, to be able to first 'run all the models' (automating that is another question for later) and then 'tidy all the models', the latter presumably generating a list of tibbles.
I've tried the following, but it throws: Error: No tidy method recognized for this list.
clean_model <- function(polr_results) {
lapply(polr_results,
broom::tidy(polr_results, conf.int = TRUE, exponentiate = TRUE) %>%
filter(coef.type=="coefficient") %>%
dplyr::arrange(-str_detect(term, 'd2sd'))
)
}
mtcars_m1 <- mtcars %>% polr(as.factor(cyl) ~ hp , data = ., Hess = TRUE)
mtcars_m2 <- mtcars %>% polr(as.factor(cyl) ~ hp + qsec , data = ., Hess = TRUE)
clean_model(c(mtcars_m1, mtcars_m2))
Upvotes: 0
Views: 500
Reputation: 11663
Another way you can do this is with purrr, putting all the different formulae you want into in a dataframe list column:
library(MASS)
library(tidyverse)
library(broom)
formula_dfs <- tibble(formula_id = 1:2,
formula = c(as.formula(as.factor(cyl) ~ hp),
as.formula(as.factor(cyl) ~ hp + qsec)))
formula_dfs
#> # A tibble: 2 x 2
#> formula_id formula
#> <int> <list>
#> 1 1 <formula>
#> 2 2 <formula>
formula_dfs %>%
mutate(polr_fit = map(formula, polr, data = mtcars, Hess = TRUE),
polr_coef = map(polr_fit, tidy, conf.int = TRUE, exponentiate = TRUE)) %>%
unnest(polr_coef) %>%
filter(coef.type=="coefficient")
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
#> # A tibble: 3 x 10
#> formula_id formula polr_fit term estimate std.error statistic conf.low
#> <int> <list> <list> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 1 <formula> <polr> hp 1.12 0.0399 2.90 1.06
#> 2 2 <formula> <polr> hp 1.13 0.0452 2.72 1.06
#> 3 2 <formula> <polr> qsec 1.18 0.369 0.448 0.538
#> # … with 2 more variables: conf.high <dbl>, coef.type <chr>
Created on 2021-05-24 by the reprex package (v2.0.0)
Your regular model output is still there in the polr_fit
column.
Upvotes: 1
Reputation: 389325
Something like this ?
library(broom)
library(tidyverse)
clean_model <- function(polr_results) {
lapply(polr_results, function(x) {
broom::tidy(x, conf.int = TRUE, exponentiate = TRUE) %>%
filter(coef.type=="coefficient")
})
}
clean_model(list(mtcars_m1, mtcars_m2))
#[[1]]
# A tibble: 1 x 7
# term estimate std.error statistic conf.low conf.high coef.type
# <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
#1 hp 1.12 0.0399 2.90 1.06 1.26 coefficient
#[[2]]
# A tibble: 2 x 7
# term estimate std.error statistic conf.low conf.high coef.type
# <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
#1 hp 1.13 0.0452 2.72 1.06 1.29 coefficient
#2 qsec 1.18 0.369 0.448 0.538 2.51 coefficient
Upvotes: 1