How to analyze the evolution of a multilevel variable with an interaction? I'm doing an analysis in a frequentist approach where the dependant variable X has 4 levels (say A, B, C, and D).
This variable was studied in 4 timeframes: 2 years (Y1 and Y2) and 2 calendar periods (P1 and P2). Here is the repartition of X depending of Y and P:

My alternative hypothesis is that this repartition is different during Y2-P2 only, independent of the year effect and the period effect. For instance, when X==C, there was an increase from P1 to P2 during Y1 but a decrease during Y2. On the other hand, when X==D, the decrease from P1 to P2 seems similar during Y1 and Y2.
Therefore, I want to fit a model as X ~ Y + P + Y:P and interpret the interaction term only. In practice, I was planning to perform a loglikelihood ratio test between this model and the nested one without this interaction term, so I can test if the interaction is significant in explaining the repartition of X.
However, as stated in my other question (What is the difference between fitting multinomal logistic regression and fitting multiple logistic regressions?), one could think of two ways to describe this:

*

*a multinomial logistic model, where I understand that an OR would represent (given all flaws in this representation) the relative risk of having a specific level instead of the reference level.

*4 logistic models on dummy binary variables, where I understand that the 4 ORs would represent the relative risk of having a specific level instead of not.

In my case as in many others, there is no reference value and expressing the risk of having D instead of A does not make any sense. I would like my conclusions to be something like: "there were a (likely) significant negative interaction in C and a non-significant interaction in B and D (and maybe A)".
Of course, one is not supposed to aim for a specific conclusion, this is just for the example.
The comments and answers to my other question showed that 4 logistic models would misestimate the probabilities. Therefore, is there a way I can express my results in the desired way using the multinomial model?
NB: I'm using nnet::multinom function in R to fit the multinomial model, a solution that can work with this or at least with R would be more than welcome!
 A: Nice question - however, I would suggest swapping your labelling so that your dependent variable is called Y and your independent variable is called X, in line with usual notation convention in statistical modelling.  It seems odd to call an independent variable Y and a dependent variable X (though I understand in your case Y stands for Year).
It's not entirely clear to me from your post what you are really interested in in terms of testing hypotheses. But it seems that your problem would involve trying to compare probabilities that the dependent variable takes a certain value.
For example, if:
Prob_A(Year1, Period1) 
Prob_A(Year1, Period2) 
Prob_A(Year2, Period1) 
Prob_A(Year2, Period2) 

are the probabilities that the dependent variable takes the value A for specific combinations of values of Year and Period, you could test hypotheses of the form:
Ho: Prob_A(Year2, Period2) - Prob_A(Year2, Period1) = Prob_A(Year1, Period2) - Prob_A(Year1, Period1)  versus
Ha: Prob_A(Year2, Period2) - Prob_A(Year2, Period1) $\neq$ Prob_A(Year1, Period2) - Prob_A(Year1, Period1)
These hypotheses would enable you to investigate whether there is a between-period difference across years with respect to the true probability that the dependent variable takes the value A.
Note that these hypotheses could be re-expressed as:
Ho: Prob_A(Year2, Period2) - Prob_A(Year2, Period1) -  (Prob_A(Year1, Period2) - Prob_A(Year1, Period1)) = 0  versus
Ha: Prob_A(Year2, Period2) - Prob_A(Year2, Period1) - (Prob_A(Year1, Period2) - Prob_A(Year1, Period1)) $\neq$ 0
Similar hypotheses can be formulated for the other levels of your dependent variable (i.e., B, C and D) and you would test all the hypotheses simultaneously by using bootstrap to build confidence intervals for quantities such as Prob_A(Year2, Period2) - Prob_A(Year2, Period1) - (Prob_A(Year1, Period2) - Prob_A(Year1, Period1)), etc., and looking to see whether these intervals exclude the value 0.
In your bootstrap procedure, you would estimate the probabilities of interest from a multinomial logistic regression model relating the dependent variable to year, period and their interaction (if the interaction is supported by the model).
A more flexible R function for fitting multinomial logistic regression models for a categorical dependent variable with unordered categories is the vglm function in the VGAM package. For example, if your model includes no interaction between Period and Year, you can fit it like this:
## Load VGAM
library(VGAM)

## Use the first level of Dependent variable 
## as the reference level for all comparisons 

model.nonparallel <- vglm(Dependent ~ Period + Year, data = yourdata,
                     family = multinomial(parallel = FALSE, 
                                          refLevel = 1))
summary(model.nonparallel)

This fits the so-called non-parallel baseline category multinomial logit model.
This vglm model assumes that, after adjusting for Year, Period has a different effect on the (i) log odds that Dependent = B rather than A, (ii) log odds that Dependent = C rather than A and (iii) log odds that Dependent = D rather than A. Similarly, after adjusting for Period, the model assumes that Year has a different effect on these odds. The non-parallel assumption refers to a situation where we allow the effect of an independent variable to vary across the log-odds correaponding to the comparisons on interest involving the levels of the dependent variables: B versus A, C versus A and D versus A.
As explained for instance in this 2017 article on On the “Poisson Trick” and its Extensions for Fitting Multinomial Regression Models by Lee et al.,
(https://arxiv.org/pdf/1707.08538.pdf),
it is possible to relax the non-parallel assumption in a non-parallel model so that we consider that some or all independent variables included in the model have effects that do NOT differ across the comparisons of interest (i.e., effects that are the same across all comparisons of interest).
A parallel baseline category multinomial logit model would be one where the effects of all independent variables included in the model are assumed to NOT vary across comparisons of interest:
model.parallel <- vglm(Dependent ~ Period + Year, data = yourdata,
                     family = multinomial(parallel = TRUE, 
                                          refLevel = 1))
summary(model.parallel)

A partial baseline category multinomial logit model would be one where the effects of some (but not all) of the independent variables included in the model are assumed to NOT vary across comparisons of interest. For example, if the effects of Period only are assumed to NOT vary across the comparisons of interest, this can be specified as:
model.partial <- vglm(Dependent ~ Period + Year, data = yourdata,
                     family = multinomial(parallel = ~ -1 + Period, 
                                          refLevel = 1))
summary(model.partial)

See this excellent post on a Unified Approach to Ordinal (cumulative) and Polytomous (multinomial) Logistic Regressions using VGAM::vglm for more vglm details: http://rstudio-pubs-static.s3.amazonaws.com/5529_7944ecfa034f4c52b8645707e48f8c6d.html.
How you conduct the bootstrap will depend on your study design. You will likely need to adjust the confidence level of your bootstrap confidence intervals for multiplicity.
If you have clear hypotheses in mind that you would like to test, such as the ones I listed here, you might want to consider keeping your interaction between Period and Year in the model and letting go of performing a test of significance of interaction.
Things get a bit more complicated if you need to compare things across levels of the dependent variable, not just within a level.
A: 
I would like my conclusions to be something like: "there were a (likely) significant negative interaction in C and a non-significant interaction in B and D (and maybe A)".

It isn't possible to test all 4 out of 4 outcome groups independently when you are modeling probabilities. Once you have determined probabilities for outcomes A, B, and C, the probability of D is completely determined. So you can examine an interaction between predictors in terms of the multinomial outcomes, but you can't assign those "significant interactions" independently among all 4 outcome groups.
That's the advantage of Isabella Ghement's suggestion (+1) to start with a proper multinomial model and then test specific (ideally, pre-specified) hypotheses about probabilities of class membership as a function of the combination of predictors of interest. Whether that's best done with the VGAM package she recommends or the nnet::multinom function you use is outside my expertise. The emmeans package in R allows expression of multinom output directly as probabilities and error estimates for all outcome classes; I don't know if it works with VGAM.
If you think it's too hard to explain to your audience the non-independence among outcome groups that's inherent in multinomial modeling, you could consider a different but related approach, as your predictors and outcome are all categorical. Log-linear analysis treats this type of problem as a multi-dimensional contingency table.
You start by modeling the numbers of observations within each of your 4 outcomes x 2 years x 2 periods = 16 combinations of outcome groups and predictor values, putting aside at first that you think of one of those dimensions as the "outcome." You model the number in each cell of the contingency table as a function of the outcome groups and your predictor values, with an implicit log link.
You start with the completely saturated model (including all interactions among outcomes and predictors) and then see which terms can be removed from the model without significantly reducing the quality of the fit. Once you have the final model, you can re-interpret in terms of your outcomes and predictor groups. This tutorial shows how to use either Poisson regression (glm) or the loglm function in the R MASS package to model an outcome as a function of categorical predictors. You might find that approach easier to explain.
