I am trying to replicate Stata's marginal effects from multinomial logit models in R but with no success. For the multinomial logit model, I used the multinom() function from the nnet package and for the marginal effects I used the margins package but the marginal_effects function seems to only display effects of a single variable. What if I want to have the marginal effects of the variable conditioned on another variable? Here is the output from Stata:
. margins, dydx(male) at(site=(1 2 3)) #male conditioned on site
Average marginal effects Number of obs = 615
Model VCE : OIM
dy/dx w.r.t. : 1.male
1._predict : Pr(insure==Indemnity), predict(pr outcome(1))
2._predict : Pr(insure==Prepaid), predict(pr outcome(2))
3._predict : Pr(insure==Uninsure), predict(pr outcome(3))
1._at : site = 1
2._at : site = 2
3._at : site = 3
------------------------------------------------------------------------------
| Delta-method
| dy/dx Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
1.male |
_predict#_at |
1 1 | -.1492951 .0728108 -2.05 0.040 -.2920016 -.0065885
1 2 | -.159346 .0723512 -2.20 0.028 -.3011517 -.0175403
1 3 | -.055138 .0875712 -0.63 0.529 -.2267745 .1164984
2 1 | .0763095 .0765406 1.00 0.319 -.0737074 .2263264
2 2 | .1747759 .0730055 2.39 0.017 .0316877 .3178641
2 3 | .0861997 .0843816 1.02 0.307 -.0791852 .2515846
3 1 | .0729855 .0516839 1.41 0.158 -.0283131 .1742842
3 2 | -.0154299 .0104982 -1.47 0.142 -.036006 .0051462
3 3 | -.0310617 .0495625 -0.63 0.531 -.1282025 .0660791
------------------------------------------------------------------------------
Note: dy/dx for factor levels is the discrete change from the base level.
My attempt to calculate the marginal effects of male using the marginal_effects function:
library(nnet)
sysdsn1$insure <- as.factor(sysdsn1$insure)
sysdsn1$male <- as.factor(sysdsn1$male)
sysdsn1$site <- as.factor(sysdsn1$site)
sysdsn1$nonwhite <- as.factor(sysdsn1$nonwhite)
sysdsn1$insure <- relevel(sysdsn1$insure, ref = "3") #set the reference level
mn0 <- multinom(insure ~ age + male*site + nonwhite, data = sysdsn1) #multinomial logit model
head(marginal_effects(mn0, variables = "male")) #this only calculate marginal effects of male, how to condition on site?
dydx_male1
1 -0.01310874
2 -0.01744213
3 0.07911846
4 -0.03386199
5 -0.01728126
6 -0.01638176
Data
Data can be downloaded from http://www.stata-press.com/data/r13/sysdsn1.dta and imported into R
