Multinomial Models for Ordinal Responses

The outcome of a response variable might be one of a restricted set of possible values. If there are only two possible outcomes, such as male and female for gender, these responses are called binary responses. If there are multiple outcomes, then they are called polytomous responses. Some examples of polytomous responses include levels of a disease (mild, medium, severe), preferred districts to live in a city, the species for a certain flower type, and so on. Sometimes there might be a natural order among the response categories. These responses are called ordinal responses.

The ordering might be inherent in the category choices, such as an individual being not satisfied, satisfied, or very satisfied with an online customer service. The ordering might also be introduced by categorization of a latent (continuous) variable, such as in the case of an individual being in the low risk, medium risk, or high risk group for developing a certain disease, based on a quantitative medical measure such as blood pressure.

You can specify a multinomial regression model that uses the natural ordering among the response categories. This ordinal model describes the relationship between the cumulative probabilities of the categories and predictor variables.

Different link functions can describe this relationship with logit and probit being the most used.

• Logit: The default link function mnrfit uses for ordinal categories is the logit link function. This models the log cumulative odds. The 'link','logit' name-value pair specifies this in mnrfit. Log cumulative odds is the logarithm of the ratio of the probability that a response belongs to a category with a value less than or equal to category j, P(y ≤ c_j), and the probability that a response belongs to a category with a value greater than category j, P(y >c_j).

  Ordinal models are usually based on the assumption that the effects of predictor variables are the same for all categories on the logarithmic scale. That is, the model has different intercepts but common slopes (coefficients) among categories. This model is called parallel regression or the proportional odds model. It is the default for ordinal responses, and the 'interactions','off' name-value pair specifies this model in mnrfit.

  The proportional odds model is

  $$\begin{array}{l}\ln \left(\frac{P(y\le c_1)}{P(y>c_1)}\right)=\ln \left(\frac{{\pi }_1}{{\pi }_2+\cdots +{\pi }_k}\right)={\alpha }_1+{\beta }_1{X}_1+{\beta }_2{X}_2+\cdots +{\beta }_p{X}_p,\\ \ln \left(\frac{P(y\le c_2)}{P(y>c_2)}\right)=\ln \left(\frac{{\pi }_1+{\pi }_2}{{\pi }_3+\cdots +{\pi }_k}\right)={\alpha }_2+{\beta }_1{X}_1+{\beta }_2{X}_2+\cdots +{\beta }_p{X}_p,\\ ⋮\\ \ln \left(\frac{P(y\le c_{k-1})}{P(y>c_{k-1})}\right)=\ln \left(\frac{{\pi }_1+{\pi }_2+\cdots +{\pi }_{k-1}}{{\pi }_k}\right)={\alpha }_{k-1}+{\beta }_1{X}_1+{\beta }_2{X}_2+\cdots +{\beta }_p{X}_p,\end{array}$$

  where π_j, j = 1, 2, ..., k, are the category probabilities.

  For example, for a response variable with three categories, there are 3 – 1 = 2 equations as follows:

  $$\begin{array}{l}\ln \left(\frac{\pi {}_1}{\pi {}_2+\pi {}_3}\right)={\alpha }_1+{\beta }_1{X}_1+{\beta }_2{X}_2+\cdots +{\beta }_p{X}_p,\\ \ln \left(\frac{\pi {}_1+\pi {}_2}{\pi {}_3}\right)={\alpha }_2+{\beta }_1{X}_1+{\beta }_2{X}_2+\cdots +{\beta }_p{X}_p.\end{array}$$

  Under the proportional odds assumption, the partial effect of a predictor variable X is invariant to the choice of the response variable category, j. For example, if there are three categories, then the coefficients express the impact of a predictor variable on the relative risk or log odds of the response value being in category 1 versus categories 2 or 3, or in category 1 or 2 versus category 3.

  Thus, a unit change in variable X₂ would mean a change in the cumulative odds of the response value being in category 1 versus categories 2 or 3, or category 1 or 2 versus category 3 by a factor of exp(β₂), given all else equal.

  You can alternatively fit a model with different intercept and slopes among the categories by using the 'interactions','on' name-value pair argument. However, using this option for ordinal models when the equal slopes model is true causes a loss of efficiency (you lose the advantage of estimating fewer parameters).

• Probit: The 'link','probit' name-value pair argument uses the probit link function which is based on a normally distributed latent variable assumption. For ordinal response variables this is also called an ordered probit model. Consider the regression model that describes the relationship of a latent variable y* of an ordinal process and a vector of predictor variables, X,

  $$y^{*}=\beta X+\epsilon ,$$

  where the error term ε has a standard normal distribution. Suppose there is the following relationship between the latent variable y* and the observed variable y:

  $$\begin{array}{l}y=c_{1}if{\alpha }_0<y^{*}\le {\alpha }_1,\\ y=c_{2}if{\alpha }_1<y^{*}\le {\alpha }_2,\\ ⋮⋮\\ y=c_{k}if{\alpha }_{k-1}<y^{*}\le {\alpha }_k,\end{array}$$

  where α₀ = – ∞ and α_k = ∞. Then, the cumulative probability of y being in category j or one of earlier categories, P(y ≤ c_j), is equal to

  $$P\left(y\le c_j\right)=P\left(y^{*}<{\alpha }_j\right)=P\left(\beta X+\epsilon <{\alpha }_j\right)=P\left(\epsilon <{\alpha }_j-\beta X\right)=\Phi \left({\alpha }_j-\beta X\right),$$

  where Φ is standard normal cumulative distribution function. Thus,

  $${\Phi }^{-1}\left(P\left(y\le c_j\right)\right)={\alpha }_j-\beta X,$$

  where α_j corresponds to the cut points of the latent variable and the intercept in the regression model. This only holds under the assumptions of a normal latent variable and parallel regression. More generally, for a response variable with k categories and multiple predictors, the ordered probit model is

  $$\begin{array}{l}{\Phi }^{-1}\left(P\left(y\le c_1\right)\right)={\alpha }_1+{\beta }_1{X}_1+\cdots +{\beta }_p{X}_p,\\ {\Phi }^{-1}\left(P\left(y\le c_2\right)\right)={\alpha }_2+{\beta }_1{X}_1+\cdots +{\beta }_p{X}_p,\\ ⋮⋮\\ {\Phi }^{-1}\left(P\left(y\le c_{k-1}\right)\right)={\alpha }_{k-1}+{\beta }_1{X}_1+\cdots +{\beta }_p{X}_p,\end{array}$$

  where P(y ≤ c_j) = π₁ + π₂ + ... + π_j.

  The coefficients indicate the impact of a unit change in the predictor variable on the likelihood of a state. A positive coefficient, β₁, for example, indicates an increase in the underlying latent variable with an increase in the corresponding predictor variable, X₁. Hence, it causes a decrease in P(y ≤ c₁) and an increase in P(y ≤ c_k).

After estimating the model coefficients using mnrfit, you can estimate the cumulative probabilities or the cumulative number in each category using mnrval with the 'type','cumulative' name-value pair option. mnrval accepts the coefficient estimates and the model statistics mnrfit returns, and estimates the categorical probabilities or the number in each category and their confidence intervals. You can specify which category or conditional probabilities or numbers to estimate by changing the value of the 'type' name-value pair argument.

References

[1] McCullagh, P., and J. A. Nelder. Generalized Linear Models. New York: Chapman & Hall, 1990.

[2] Long, J. S. Regression Models for Categorical and Limited Dependent Variables. Sage Publications, 1997.

[3] Dobson, A. J., and A. G. Barnett. An Introduction to Generalized Linear Models. Chapman and Hall/CRC. Taylor & Francis Group, 2008.