1

I am trying to display both unstandardized and standardized coefficients in a table using modelsummary (incredibly useful package, by the way) in a wide format. Calculating the standardized coefficients works fine and they are displayed in the table. However, I would like to also display the unstandardized coefficients with their respective SE as well. My code so far looks something like this:

models = list(
  "Model X" <- lm(mpg ~ hp + factor(cyl), data = mtcars),
  "Model Y" <- lm(mpg ~ hp + factor(cyl) + qsec, data = mtcars))
names(models) <-  c("Model X", "Model Y")
test <- modelsummary(models, 
                      standardize = "basic",
                      shape = term~model + statistic,
                      statistic = c("statistic", "p.value"),
                      estimate = "{estimate} ({std.error})",
                      stars = TRUE,
                      fmt = 2,
                      col.names = c("", "Std. Est.  (SE)", "t", "p", "Std. Est (SE)", "t", "p"),
                      gof_omit = 'AIC|BIC|Log.Lik.|RMSE',
                      notes = "Note: Model Y includes control variables"
                      
) %>%kable_styling(font_size = 8) %>% row_spec(0, italic = T)
test

At the moment, the resulting table only shows the standardized coefficients. In addition, two other questions I had:

  • Stars are currently not showing, is there any possibility to show them after the p-values?
  • Is it possible to give three decimals for p-values, but keep two decimals for all other columns?
statleo
  • 81
  • 1
  • 7

2 Answers2

4

As noted in the documentation ?modelsummary, additional arguments not explicitly supported by the modelsummary() function are pushed forward to the parameters::parameters() function. standardize is one of those arguments, and it will be applied to all models.

One workaround is to extract the models into an intermediate representation by setting output="modelsummary_list". Then, we can feed those back to modelsummary to create the final table:

library(modelsummary)
mod <- lm(mpg ~ hp + factor(cyl), data = mtcars)

mod1 <- modelsummary(mod, output = "modelsummary_list")
mod2 <- modelsummary(mod, output = "modelsummary_list", standardize = "basic")

models <- list(
    "Unstandardized" = mod1,
    "Standardized" = mod2)

modelsummary(
    models,
    output = "markdown",
    stars = TRUE,
    shape = term + model ~ statistic,
    statistic = c("std.error", "p.value"),
    fmt = list(estimate = 3, p.value = 2))
Est. S.E. p
(Intercept) Unstandardized 28.650*** 1.588 0.00
Standardized 0.000*** 0.000 0.00
hp Unstandardized -0.024 0.015 0.13
Standardized -0.273 0.175 0.13
factor(cyl)6 Unstandardized -5.968** 1.639 0.00
Standardized -0.416** 0.114 0.00
factor(cyl)8 Unstandardized -8.521** 2.326 0.00
Standardized -0.713** 0.195 0.00

Note: ^^ + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

See the fmt argument above, and refer to the ?modelsummary documentation for formatting with different digits for different statistics.

In the table above, you’ll see that the stars appear after the estimates. You can make them appear after the p values instead by using glue strings, but there are a couple quirks that make this process suboptimal. I opened two issues on Github to track progress on improving the situation there:

modelsummary(
    models,
    shape = term + model ~ statistic,
    estimate = "{estimate}",
    statistic = "{p.value}{stars}",
    fmt = list(estimate = 3, p.value = 2))
Est. {p.value}{stars}
(Intercept) Unstandardized 28.650 0.00***
Standardized 0.000 0.00***
hp Unstandardized -0.024 0.13
Standardized -0.273 0.13
factor(cyl)6 Unstandardized -5.968 0.00**
Standardized -0.416 0.00**
factor(cyl)8 Unstandardized -8.521 0.00**
Standardized -0.713 0.00**
Vincent
  • 15,809
  • 7
  • 37
  • 39
  • Thank you very much for your fast reply and the clarification. Your workaround looks quite good. However, I was hoping to "simply" add an additional column with the (un-)standardized coefficients. Do you have any idea of how to achieve this? – statleo Aug 09 '22 at 12:23
  • 1
    What do you mean by "simply", exactly? In Section 4.3 of my JSS article I show how to add Bonferroni-corrected p values by defining a `tidy_custom.lm()` method. You could do the same with standardized estimates. https://www.jstatsoft.org/article/view/v103i01 – Vincent Aug 09 '22 at 13:52
  • With "simply" I just wanted to emphasize that even though something might look simple (i.e., adding a new column to the table), it is not really simple. Thanks for linking the article, the tidy_custom.lm() helped a lot and I was able to solve my issue (also looking at your git repository). – statleo Aug 09 '22 at 18:53
  • 1
    Got it. Yeah, sometimes things get more complicated. Honestly, in a lot of cases it makes sense to just build a data frame (perhaps with the help of `get_estimates(mod, standardize="basic")`, and the feed that data frame to `add_column`. But I'm glad you found a solution! – Vincent Aug 09 '22 at 19:17
  • Sorry for annoying you again, but I tried using my solution with IV regressions (ivreg) and now I get the error that object 'estimate.unst' was not found. If I follow the steps of the tidy_custom.lm function and enter the IV model, I get a df with the respective columns that looks just like the one for a regular lm. Any ideas on what could be the issue? – statleo Aug 10 '22 at 12:05
  • 1
    Did you redefine a new `tidy_custom.ivreg()` method? Different class; different method. Also, if you execute the same code line by line, can you get the `s[, 1]` in the same way as in the previous kind of model? – Vincent Aug 10 '22 at 13:42
  • 1
    Worked perfectly, Haven't worked with these methods and functions before, but they are insanely helpful. Thank you so much! – statleo Aug 10 '22 at 14:10
0

I solved the issue by using this tidy_custom.lm solution as suggested by Vincent:

tidy_custom.lm <- function(x, ...) {
  s <- summary(x)$coefficients
  out <- data.frame(
    term = row.names(s),
    estimate.unstd = s[,1],
    se.unstd = s[,2])
  return(out)
}

In the modelsummary, I adjusted the following commands:

statistic = c("{estimate} ({std.error})", "statistic", "p.value"),
estimate = "{estimate.unstd} ({se.unstd}) {stars}"
statleo
  • 81
  • 1
  • 7