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

Question

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:

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.

score 2 · Accepted Answer · answered Jul 13 '22 at 17:33

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:

Display term names and contrast identifiers in separate columns, and rename them using the group_map argument.
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

Thank you @Vincent, the merging of the two columns indeed works well for me. I discovered that the reason why I could produce the table in the image I included was because I had dichotomized my variables completely with 0-1 (as suggested by @FilippeO.Gonçalves). The new features certainly make this process more efficient. Thank you. — Daniel Sánchez, Jul 16 '22 at 05:25

score 1 · Answer 2 · answered Jun 30 '22 at 00:17

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

Thank you. It actually is not necessary for me to present those stats, so that isn't a problem. I understood the thing about the column, I just wish I could edit the stuff in the column, changing dy/dx to something else, for instance. I'm also surprised that *marginaleffects* objects are not supported as they are mantained by Vincent too, and there is an example using that packge in the *modelsummary* vignettes. — Daniel Sánchez, Jun 30 '22 at 05:11
Glad I've helped! Yeah, that incompatibility surprised me too. — Filippe O. Gonçalves, Jun 30 '22 at 15:00
I agree that the interop between the two package is not optimal. The warnings are a known problem, and the solution should be coming soon. Combining the two columns more easily would be a nice feature. I'm not exactly sure what's the best user interface for this, but I'll think about it. I'm on vacation now, but will look into this when I return. Thanks for using the packages! — Vincent, Jun 30 '22 at 20:12
@DanielSánchez I have added new functionality to the development version of the package and added an answer to illustrate the new functions. If you get a chance to try it, let me know if it works for you. — Vincent, Jul 13 '22 at 17:34
@FilippeO.Gonçalves the new answer may also be of interest to you. — Vincent, Jul 13 '22 at 17:35

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

2 Answers2