R lsmeans adjust multinomial regression model with interaction term

Question

I want to run a multinomial model with mblogit(), which the dependent variable (Treatment) is a three level categorical variable ("treatment1"(reference), "treatment2" and "treatment3"), and the predictors (Exposure) include a 4 level categorical variable("low"(reference), "medium", "high", "very high") and Year information (2016(reference), 2017, 2018, 2019,2020). The model also includes a random effect (Hospital), see below:

model <- mblogit(Treatment ~ Exposure + Year + Exposure*Year, 
                 data = my_df, catCov ="free", random =  ~1|Hospital)

Now I want to get the lsmeans for each Exposure level by Years like below:

Year = 2016: Exposure Estimate SE P-value
low -- -- --- medium -- --- --- high -- --- --- very high -- ---- ----

Year = 2017: Exposure Estimate SE P-value
low -- -- --- medium -- --- --- high -- --- --- very high -- ---- ----

I tried this code to get lsmeans():

lsmeans(model, "Exposure", by = "Year", adjust = "tukey")

and I got these results which is not what I expect:

Year = 2016:
 Exposure    prob       SE      df    asymp.LCL  asymp.UCL
 low         0.333      NaN     Inf       NaN       NaN
 medium      0.333   0.00e+00   Inf     0.333     0.333
 high        0.333      NaN     Inf       NaN       NaN
 very high   0.333   3.91e-10   Inf     0.333     0.333

Year = 2017:
 Exposure    prob       SE      df    asymp.LCL  asymp.UCL
 low         0.333      NaN     Inf       NaN       NaN
 medium      0.333   0.00e+00   Inf     0.333     0.333
 high        0.333      NaN     Inf       NaN       NaN
 very high   0.333   3.91e-10   Inf     0.333     0.333

Year = 2018