I'd like to be able to use the LRstats function from the vcdExtra package to compare 50 models. However, this function appears to be set up such that it needs each model object provided individually. If I have a list of 50 models, how can I do this without having to specify each model manually i.e. test[[1]], test[[2]], test[[3]]...etc?
library(vcdExtra)
library(dplyr)
test <- data.frame(
y = rbinom(2500, 1, 0.5),
x1 = rnorm(2500, mean = 0, sd = 1),
x2 = rnorm(2500, mean = 1, sd = 3),
z = rep(seq(from = 400, length.out = 50, by = 400), times = 50))
test_list <- group_split(test, z)
names(test_list) <- unique(test$z)
test_models <- lapply(test_list, function(x) glm(y ~ x1 x2, data = x, family = "binomial"))
test_models2 <- glmlist(test_models)
>Warning message:
In glmlist(test_models) :
Objects test_models removed because they are not glm objects
test_LR <- LRstats(test_models2)
>NULL
CodePudding user response:
The problem was in the way function glmlist
was used. glmlist
expects one or more objects inheriting from class "glm"
or "loglm"
. From the documentation:
Any objects which do not inherit the appropriate class glm or loglm are excluded, with a warning.
So call the function on each list member with do.call
.
library(vcdExtra)
#> Loading required package: vcd
#> Loading required package: grid
#> Loading required package: gnm
test <- data.frame(
y = rbinom(2500, 1, 0.5),
x1 = rnorm(2500, mean = 0, sd = 1),
x2 = rnorm(2500, mean = 1, sd = 3),
z = rep(seq(from = 400, length.out = 50, by = 400), times = 50))
test_list <- dplyr::group_split(test, z)
names(test_list) <- unique(test$z)
test_models <- lapply(test_list, function(x) glm(y ~ x1 x2, data = x, family = "binomial"))
test_models2 <- do.call(glmlist, test_models)
test_LR <- LRstats(test_models2)
head(test_LR)
#> Likelihood summary table:
#> AIC BIC LR Chisq Df Pr(>Chisq)
#> 400 72.136 77.872 66.136 47 0.03420 *
#> 800 68.608 74.344 62.608 47 0.06340 .
#> 1200 74.842 80.578 68.842 47 0.02056 *
#> 1600 73.545 79.281 67.545 47 0.02634 *
#> 2000 73.514 79.250 67.514 47 0.02649 *
#> 2400 70.524 76.260 64.524 47 0.04564 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Created on 2022-03-31 by the reprex package (v2.0.1)