1

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:

barplot

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!

kjetil b halvorsen
  • 63,378
  • 26
  • 142
  • 467
Dan Chaltiel
  • 1,089
  • 12
  • 25

2 Answers2

2

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.

Isabella Ghement
  • 18,164
  • 2
  • 22
  • 46
  • 1
    Thank you very much for your answer. Maybe I was not clear enough. I plan to perform a loglikelihood ratio test between my model and the nested one without the interaction term. This way, I should be able to test if the interaction is significant in explaining the repartition of `X`. (Yes, it bothers me as well to call it `X` but the year was `Y` and I'd have to make the plot again for which I didn't keep the code, sorry for that!). – Dan Chaltiel Nov 16 '20 at 08:21
  • 1
    The bootstrap procedure is very clever, I'll definitely try to give it a try! However, I'm not sure to understand, what would `glmm` bring to `nnet::multinom`? Also, are you sure that the proportional odds assumption is needed for multinomial models? – Dan Chaltiel Nov 16 '20 at 08:24
  • 1
    Keep all code, Dan - better safe than sorry! The VGAM package I recommended after editing my answer does have an lrtest() function for performing a likelihood ratii test involving two competing, nested models. However, if you have clear hypotheses you are interested in testing, my feeling is that you should keep the interaction term in the model and proceed to test those hypotheses. Planned comparisons typically justify skipping an (omnibus) test of significance of interaction. – Isabella Ghement Nov 16 '20 at 20:01
  • I cleared things up in my answer - the use of proportional odds was not warranted; what I wanted to say was “non-parallel” versus “parallel”/“partial” assumption, but I answered your question late in the evening and got muddled in the terminology. – Isabella Ghement Nov 16 '20 at 20:04
  • 1
    This is an incredible piece of an answer, thanks a lot! Actually, I thought about the LRT to test whether the interaction term was significant, I feel like this is my very hypothesis. These `vglm` models are also very interesting for many hypothesis. However, I feel like they unfortunately don't adress my problem as, just like `nnet::multinom`, they need a reference level which does not make sense in my case. – Dan Chaltiel Nov 17 '20 at 11:04
  • You’re welcome, Dan! These multinomial logistic regression models can be interpreted on 3 different scales - log odds, odds (via odds ratios) and probabilities. In your initial post, you mentioned that “your alternative hypothesis is that this repartition is fifferent, etc.” That doesn’t elucidate what scale you are actually interested in when formulating your hupotheses. If you are interested in the probability scale, you would formulate your hypotheses in terms of probabilities, much like I did in my answer. – Isabella Ghement Nov 17 '20 at 16:15
  • With hypotheses formulated in terms of probabilities, it shouldn’t really matter that you fit your model by picking a reference level for your response variable (e.g., A) - that reference level would only come into play if you were interested in the log odds or odds scale. All you would be interested in is using the model to estimate the probabilities of interest from many bootstrap samples and construct bootstrap confidence intervals for the linear combinations of probsbilities you care about, with the coefficients in those combinations being 1 or -1. – Isabella Ghement Nov 17 '20 at 16:19
  • The coefficient of the interaction between Period and Year “lives” on the log odds scale - if you are interested in the probability scale, this coefficient will enter the pribabilities of interest (e.g., probability that response value is A in Period 1 of Year 1) in a nonlinear fashion. I doubt testing whether or not this coefficient is 0 on the (linear!) log odds scale will tell you anything useful about probability-based hypotheses such as the ones I listed in my answer. – Isabella Ghement Nov 17 '20 at 16:23
  • I meant *anything* **directly** *useful*... – Isabella Ghement Nov 17 '20 at 17:36
  • Let us [continue this discussion in chat](https://chat.stackexchange.com/rooms/116329/discussion-between-dan-chaltiel-and-isabella-ghement). – Dan Chaltiel Nov 17 '20 at 17:39
1

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.

EdM
  • 57,766
  • 7
  • 66
  • 187