I work mostly with binary and ordinal outcomes. So, I rarely use multinomial logistic regression but yesterday I had to, and it proved tricky to create tables for model results. Basically, multinomial logistic regression generalizes logistic regression to nominal variables with multiple categories. It is preferred when these categories are unordered but also used for ordinal outcomes when proportional odds assumption is violated. It is easy to estimate, but not as straightforward to interpret as binary logistic regression. And it might not be practical when your outcome has many categories. In R, you can use multinom function from {nnet} package to estimate multinomial logistic regression models:1
1 I use the data from the example on IDRE UCLA: https://stats.idre.ucla.edu/r/dae/multinomial-logistic-regression/.
df <-read.dta("https://stats.idre.ucla.edu/stat/data/hsbdemo.dta") multinom_fit <-multinom(prog ~ ses + female + write, data = df, trace=FALSE) # trace=FALSE to suppress convergence info.summary(multinom_fit)
There are no \(p\)-values on this output but you can calculate them manually as explained in the IDRE UCLA example (see the footnote). That is not our issue now. We want to create a table for these results in a publishable format. And that is the tricky part. {Stargazer} provides one (recommended) method to do so:
library(stargazer)stargazer(multinom_fit, type ="html")
Dependent variable:
academic
vocation
(1)
(2)
sesmiddle
0.519
0.926*
(0.451)
(0.502)
seshigh
1.147**
0.328
(0.523)
(0.663)
femalefemale
-0.080
0.526
(0.397)
(0.463)
write
0.059***
-0.066***
(0.022)
(0.025)
Constant
-2.854**
2.489**
(1.166)
(1.185)
Akaike Inf. Crit.
377.718
377.718
Note:
p<0.1; p<0.05; p<0.01
You can also export the table to \(\LaTeX\), but not Word (don’t forget to add results='asis' to code chunk if you use rmarkdown document). Next, you can use tab_model from {sjPlot}:
library(sjPlot)tab_model(multinom_fit)
prog
Predictors
Odds Ratios
CI
p
Response
(Intercept)
0.06
0.01 – 0.58
0.015
academic
ses [middle]
1.68
0.69 – 4.09
0.251
academic
ses [high]
3.15
1.12 – 8.84
0.029
academic
female [female]
0.92
0.42 – 2.02
0.841
academic
write
1.06
1.01 – 1.11
0.009
academic
(Intercept)
12.05
1.16 – 124.67
0.037
vocation
ses [middle]
2.52
0.94 – 6.79
0.066
vocation
ses [high]
1.39
0.38 – 5.13
0.621
vocation
female [female]
1.69
0.68 – 4.22
0.257
vocation
write
0.94
0.89 – 0.98
0.010
vocation
Observations
200
R2 / R2 adjusted
0.124 / 0.119
tab_model is great to create html tables. You can then drag and drop the html file to a Word document, and tweak a little to make it publication-ready. But it does not export to \(\LaTeX\). Moreover, you might want a wide table (outcome categories side-by-side) rather than a long one.
Other options are available, e.g., {texreg} package. It provides functions to export in all three formats mentioned above. But they create long tables.
This could do the job if you are okay with a long table but it requires some work. Unfortunately, I wanted a wide table with outcome categories on separate columns and I needed to export to Word.
The solution comes from the {modelsummary} package. Using the shape argument in the modelsummary function, we could create tables with outcome categories side-by-side.
library(modelsummary)library(flextable)modelsummary(multinom_fit,stars =TRUE, exponentiate =TRUE, shape = model + term ~ response, output ="flextable") %>%autofit()
academic
vocation
(1)
(Intercept)
0.058*
12.051*
(0.067)
(14.274)
sesmiddle
1.680
2.524+
(0.757)
(1.266)
seshigh
3.150*
1.388
(1.647)
(0.920)
femalefemale
0.923
1.693
(0.366)
(0.784)
write
1.061**
0.937**
(0.024)
(0.023)
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001
If you have more than one model, as usually is the case, you might want to specify shape argument in a different way:
library(modelsummary)library(flextable)multinom_fit1 <-multinom(prog ~ ses + female + write, data = df, trace=FALSE)multinom_fit2 <-multinom(prog ~ ses + female + write + math, data = df, trace=FALSE)modelsummary(list(multinom_fit1, multinom_fit2),stars =TRUE, exponentiate =TRUE, shape = term ~ model + response, output ="flextable") %>%autofit()
(1) / academic
(1) / vocation
(2) / academic
(2) / vocation
(Intercept)
0.058*
12.051*
0.012**
42.394*
(0.067)
(14.274)
(0.017)
(64.709)
sesmiddle
1.680
2.524+
1.436
2.700*
(0.757)
(1.266)
(0.669)
(1.367)
seshigh
3.150*
1.388
2.578+
1.552
(1.647)
(0.920)
(1.395)
(1.044)
femalefemale
0.923
1.693
1.167
1.480
(0.366)
(0.784)
(0.480)
(0.699)
write
1.061**
0.937**
1.016
0.952+
(0.024)
(0.023)
(0.027)
(0.028)
math
1.075**
0.959
(0.030)
(0.032)
Num.Obs.
200
200
R2
0.124
0.165
R2 Adj.
0.119
0.160
AIC
377.7
365.0
BIC
410.7
404.6
RMSE
0.42
0.41
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001
As you can see, I also specify output argument as "flextable". Thanks to {modelsummary} package, we can use {gt}, {kableExtra}, {huxtable}, and {flextable} to customize the appearance of the table. I prefer {flextable} because it makes adding the table into a Word document quite easy as explained here.
We also need to make other changes: adding title, labels, spanner, confidence intervals, goodnes-of-fit statistics, etc.2 Luckily, there are various helper functions to make these changes. I highly recommend checking the websites for {modelsummary} and {flextable}.
2 When you use odds ratios, and associated standard errors and confidence intervals, you need to be extra careful: https://github.com/tidymodels/broom/issues/422. For some, standard error of an odds ratio does not make much sense.
References
Arel-Bundock, Vincent. 2021. Modelsummary: Summary Tables and Plots for Statistical Models and Data: Beautiful, Customizable, and Publication-Ready. https://CRAN.R-project.org/package=modelsummary.
Hlavac, Marek. 2018. Stargazer: Well-Formatted Regression and Summary Statistics Tables. Bratislava, Slovakia: Central European Labour Studies Institute (CELSI). https://CRAN.R-project.org/package=stargazer.
Iannone, Richard, Joe Cheng, and Barret Schloerke. 2021. Gt: Easily Create Presentation-Ready Display Tables. https://CRAN.R-project.org/package=gt.
Leifeld, Philip. 2013. “texreg: Conversion of Statistical Model Output in R to LaTeX and HTML Tables.”Journal of Statistical Software 55 (8): 1–24. http://dx.doi.org/10.18637/jss.v055.i08.