cliu
cliu

Reputation: 965

How to add interaction terms in multinomial regression

I am using the mlogit function from the mlogit package to run a multinomial logit regression. I am not sure how to add interaction terms into my model. Here is a toy dataset and my attempt to add interactions:

library(mlogit)
data <- data.frame(y=sample(1:3, 24, replace = TRUE), 
        x1 = c(rep(1,12), rep(2,12)),
        x2 = rep(c(rep(1,4), rep(2,4), rep(3,4)),2),
        x3=rnorm(24),
        z1 = sample(1:10, 24, replace = TRUE))

m0 <- mlogit(y ~ 0|x1 + x2 + x3 + z1, shape = "wide", data = data) #model with only main effects
m1 <- mlogit(y ~ 0|(x1 + x2 + x3 + z1)^2, shape = "wide", data = data) #model assuming with all possible 2-way interactions?

The output from summary(m1) shows:

Coefficients :
               Estimate Std. Error z-value Pr(>|z|)
(Intercept):2  86.41088  164.93831  0.5239   0.6003
(Intercept):3  62.43859  163.57346  0.3817   0.7027
x1:2          -32.27065   82.62474 -0.3906   0.6961
x1:3            0.24661   84.07429  0.0029   0.9977
x2:2          -75.09247   81.36496 -0.9229   0.3561
x2:3          -85.16452   81.40983 -1.0461   0.2955
x3:2          113.11778  119.15990  0.9493   0.3425
x3:3          112.77622  117.74567  0.9578   0.3382
z1:2           11.18665   22.32508  0.5011   0.6163
z1:3           13.15552   22.26441  0.5909   0.5546
x1:2           34.01298   39.66983  0.8574   0.3912
x1:3           32.19141   39.48373  0.8153   0.4149
x1:2          -53.86747   59.75696 -0.9014   0.3674
x1:3          -47.97693   59.09055 -0.8119   0.4168
x1:2           -6.98799   11.29920 -0.6185   0.5363
x1:3          -10.41574   11.52313 -0.9039   0.3660
x2:2            0.59185    6.68807  0.0885   0.9295
x2:3            2.63458    4.94419  0.5329   0.5941
x2:2            0.80945    2.03769  0.3972   0.6912
x2:3            2.60383    2.21878  1.1735   0.2406
x3:2           -0.64112    1.64678 -0.3893   0.6970
x3:3           -2.14289    1.98436 -1.0799   0.2802

the first column is not quite clear to me what specific interactions were outputted. Any pointers will be greatly appreciated!

Upvotes: 0

Views: 1938

Answers (2)

TeeCee
TeeCee

Reputation: 15

This seems to be an issue with the way the mlogit package prints coefficient names. The number of coefficients shown is correct, but it doesn't correctly display the names of the interaction terms. The nnet::multinom() workaround is fine, but isn't sufficient if e.g., the user is looking for clustered standard errors -- the sandwich package has a vcovCL() method for mlogit models but not nnet models.

Upvotes: 1

Keith McNulty
Keith McNulty

Reputation: 972

This might be a clearer way to do it:

library(dplyr)
library(broom)
library(nnet)

multinom(formula = y ~ (x1 + x2 + x3 + z1)^2, data = data) %>% 
  tidy()

# A tibble: 22 x 6
   y.level term        estimate std.error statistic p.value
   <chr>   <chr>          <dbl>     <dbl>     <dbl>   <dbl>
 1 2       (Intercept)   -158.       247.   -0.640    0.522
 2 2       x1            -388.       247.   -1.57     0.116
 3 2       x2             -13.4      248.   -0.0543   0.957
 4 2       x3             120.       334.    0.360    0.719
 5 2       z1             173.       968.    0.179    0.858
 6 2       x1:x2          337.       248.    1.36     0.174
 7 2       x1:x3           40.2      334.    0.120    0.904
 8 2       x1:z1          -53.8      968.   -0.0555   0.956
 9 2       x2:x3         -137.      1018.   -0.135    0.893
10 2       x2:z1          -76.6      910.   -0.0841   0.933
# … with 12 more rows

Upvotes: 1

Related Questions