CODEWITHSUNDEEP
rcoefficientsmodelsummary
Daniel Sánchez
Daniel Sánchez

Reputation: 45

Getting weird output when using the modelsummary package in R to create an average partial effects table

I am trying to present an average partial effects table of multiple binary-response models. Back in December last year, I had no issue calling modelsummary in a list of objects of the class "margins" "data.frame" to produce a table like the one below: enter image description here

However, when I try to call modelsummary on a list of margins or marginaleffects objects, I get the following:

# Working example

df<-mtcars

df$cyl<-as.factor(df$cyl)

library(modelsummary)
library(tidyverse)
library(marginaleffects)

# Binary response Models

model1<-glm(am ~ mpg + cyl,
            data = df,
            family = quasibinomial(link = 'logit'))

model2<-glm(am ~ mpg + cyl,
            data = df,
            family = quasibinomial(link = 'probit'))

model3<-glm(am ~ wt + cyl,
            data = df,
            family = quasibinomial(link = 'logit'))

model4<-glm(am ~ wt + cyl,
            data = df,
            family = quasibinomial(link = 'probit'))

models<-list(model1,model2,model3,model4)

mfx<-lapply(models, marginaleffects)

modelsummary(mfx,
             output = 'markdown')
Model 1 Model 2 Model 3 Model 4
mpg 0.056 0.057
(0.027) (0.026)
cyl 0.093 0.091 0.120 0.121
0.093 0.091 0.120 0.260
0.093 0.091 0.253 0.121
0.093 0.091 0.253 0.260
0.093 0.102 0.120 0.121
0.093 0.102 0.120 0.260
0.093 0.102 0.253 0.121
0.093 0.102 0.253 0.260
0.097 0.091 0.120 0.121
0.097 0.091 0.120 0.260
0.097 0.091 0.253 0.121
0.097 0.091 0.253 0.260
0.097 0.102 0.120 0.121
0.097 0.102 0.120 0.260
0.097 0.102 0.253 0.121
0.097 0.102 0.253 0.260
(0.169) (0.174) (0.058) (0.053)
(0.169) (0.174) (0.058) (0.063)
(0.169) (0.174) (0.068) (0.053)
(0.169) (0.174) (0.068) (0.063)
(0.169) (0.236) (0.058) (0.053)
(0.169) (0.236) (0.058) (0.063)
(0.169) (0.236) (0.068) (0.053)
(0.169) (0.236) (0.068) (0.063)
(0.238) (0.174) (0.058) (0.053)
(0.238) (0.174) (0.058) (0.063)
(0.238) (0.174) (0.068) (0.053)
(0.238) (0.174) (0.068) (0.063)
(0.238) (0.236) (0.058) (0.053)
(0.238) (0.236) (0.058) (0.063)
(0.238) (0.236) (0.068) (0.053)
(0.238) (0.236) (0.068) (0.063)
wt -0.533 -0.556
(0.080) (0.072)
Num.Obs. 32 32 32 32
F 2.157 2.707 4.417 5.832
RMSE 0.39 0.39 0.27 0.27

I get the following warning messages:

--- 
Objects of class 'NULL' are currently not supported.
Objects of class 'NULL' are currently not supported.
Objects of class 'NULL' are currently not supported.
Objects of class 'NULL' are currently not supported.
Warning message:
There are duplicate term names in the table. The `shape` argument of the `modelsummary` function can be used to print related terms together, and to label them appropriately. You can find the group identifier to use in the `shape` argument by calling `get_estimates()` on one of your models. Candidate group identifiers include: type, contrast. See `?modelsummary` for details.

I tried listening to the warning message and added a shape argument as shape = ~ term + contrast ~ model, which results in a better table. But I still get the "objects of class NULL" warning and I also get another new column on my table:

Model 1 Model 2 Model 3 Model 4
mpg dY/dX 0.056 0.057
(0.027) (0.026)
cyl 6 - 4 0.097 0.091 0.120 0.121
(0.169) (0.174) (0.058) (0.053)
8 - 4 0.093 0.102 0.253 0.260
(0.238) (0.236) (0.068) (0.063)
wt dY/dX -0.533 -0.556
(0.080) (0.072)
Num.Obs. 32 32 32 32
F 2.157 2.707 4.417 5.832
RMSE 0.39 0.39 0.27 0.27

I would like to know if there is any way to get rid of that warning message, and also how to ommit/edit the new column that appears in the table with the shape argument.

Upvotes: 0

Views: 1166

Answers (2)

Vincent
Vincent

Reputation: 17725

This answer uses the development version of modelsummary, which can be installed as follows (make sure you restart R after install):

library(remotes)
install_github("vincentarelbundock/modelsummary")

There are two main approaches:

  1. Display term names and contrast identifiers in separate columns, and rename them using the group_map argument.
  2. Use the shape argument to combine the term and contrast identifiers.

First, note that contrasts and slopes are uniquely identified by two columns: term and contrast

library(modelsummary)
library(marginaleffects)

mod <- list(
    glm(am ~ factor(cyl), data = mtcars, family = binomial),
    glm(am ~ factor(cyl) + mpg, data = mtcars, family = binomial))

mfx <- lapply(mod, marginaleffects)

tidy(mfx[[1]])
#>       type term contrast   estimate std.error statistic      p.value
#> 1 response  cyl    6 - 4 -0.2987013 0.2302534 -1.297272 0.1945376094
#> 2 response  cyl    8 - 4 -0.5844156 0.1636389 -3.571373 0.0003551141
#>     conf.low  conf.high
#> 1 -0.7499897  0.1525871
#> 2 -0.9051419 -0.2636893

We can display terms and contrasts in separate columns using the shape argument:

modelsummary(
    mfx,
    shape = term + contrast ~ model
)
Model 1 Model 2
cyl 6 - 4 -0.299 0.097
(0.230) (0.166)
8 - 4 -0.584 0.093
(0.164) (0.233)
mpg dY/dX 0.056
(0.027)
Num.Obs. 32 32
AIC 39.9 37.4
BIC 44.3 43.3
Log.Lik. -16.967 -14.702
F 3.691 2.236
RMSE 0.42 0.39

We can rename and reorder the groups with the group_map argument:

modelsummary(
    mfx,
    group_map = c("8 - 4" = "High - Low", "6 - 4" = "Mid - Low", "dY/dX" = "Slope"),
    shape = term + contrast ~ model
)
Model 1 Model 2
cyl High - Low -0.584 0.093
(0.164) (0.233)
Mid - Low -0.299 0.097
(0.230) (0.166)
mpg Slope 0.056
(0.027)
Num.Obs. 32 32
AIC 39.9 37.4
BIC 44.3 43.3
Log.Lik. -16.967 -14.702
F 3.691 2.236
RMSE 0.42 0.39

We can display terms and group identifiers in the same column by including an interaction in the shape formula:


modelsummary(
    mfx,
    shape = term : contrast ~ model
)
Model 1 Model 2
cyl 6 - 4 -0.299 0.097
(0.230) (0.166)
cyl 8 - 4 -0.584 0.093
(0.164) (0.233)
mpg dY/dX 0.056
(0.027)
Num.Obs. 32 32
AIC 39.9 37.4
BIC 44.3 43.3
Log.Lik. -16.967 -14.702
F 3.691 2.236
RMSE 0.42 0.39

We can then rename or omit using the coef_* arguments:

modelsummary(
    mfx,
    shape = term : contrast ~ model,
    coef_rename = function(x) gsub("dY/dX", "(Slope)", x)
)
Model 1 Model 2
cyl 6 - 4 -0.299 0.097
(0.230) (0.166)
cyl 8 - 4 -0.584 0.093
(0.164) (0.233)
mpg (Slope) 0.056
(0.027)
Num.Obs. 32 32
AIC 39.9 37.4
BIC 44.3 43.3
Log.Lik. -16.967 -14.702
F 3.691 2.236
RMSE 0.42 0.39

Upvotes: 2

Filippe O. Gon&#231;alves
Filippe O. Gon&#231;alves

Reputation: 327

There's a lot to discuss here:

1. The new column that appeared was necessary because of the variable cyl. Since the glm function dummifies categorical variables internally, you end up with 2 variables out of cyl (it have 3 categories, so 3 minus 1 to avoid perfect multicolinearity). The table have to distinguish them in some way.

2. The default style of summaries does not yet support a "marginaleffects" model class. It's a potential issue that you can open here, helping the package development.

3. Unfortunately, I couldn't figured out how to include the F-statistic and the RMSE metric to the final table, so probably my answer is not complete in terms of desired output. This change of function behavior was weird.

# compatible style
options(modelsummary_get = "broom")

df <- mtcars

# manually dummifying
df$cyl4 <- ifelse(df$cyl == 4, 1, 0)
df$cyl6 <- ifelse(df$cyl == 6, 1, 0)

# (just avoiding excessive typing)
df <- subset(df, select = c("am", "cyl4", "cyl6", "mpg", "wt"))
glm2 <- function(x, y) glm(x, data = df, family = quasibinomial(y))

model1 <- glm2(am ~ . - wt,  "logit")
model2 <- glm2(am ~ . - wt,  "probit")
model3 <- glm2(am ~ . - mpg, "logit")
model4 <- glm2(am ~ . - mpg, "probit")

models <- list(model1, model2, model3, model4)
mfx <- lapply(models, marginaleffects::marginaleffects)

modelsummary::modelsummary(mfx, output = "markdown", gof_map = "nobs")
Model 1 Model 2 Model 3 Model 4
cyl4 -0.106 -0.116 -0.365 -0.389
(0.300) (0.299) (0.137) (0.132)
cyl6 0.005 -0.011 -0.154 -0.164
(0.205) (0.202) (0.100) (0.093)
mpg 0.056 0.057
(0.027) (0.026)
wt -0.533 -0.556
(0.080) (0.072)
Num.Obs. 32 32 32 32

Upvotes: 1

Related Questions