0

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!

cliu
  • 933
  • 6
  • 13

1 Answers1

1

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
Keith McNulty
  • 952
  • 1
  • 6
  • 17