6.6 Ordinal response variable
When the response variable still has a small number of classes but exhibits some order between those classes is called an ordinal response variable. For instance, the income with the levels, low, medium, and high, can be considered as an ordinal one. To predict such a variable, one can use the multinomial model, discussed earlier, by ignoring the ordering property. However, since this property can be critical for most situations, predictions will be less powerful. Researchers proposed a bunch of models for this type of model. In this subsection, we will discuss the most used ones.
6.6.1 Cumulative logit model
It was introduced first by (Walker.SH and Duncan.DB 1967). The idea behind this model is to convert the problem into response binary models (described in 6.4) where their number is equal to the number of the response classes.
For simplicity, suppose we have data with an ordered outcome \(y\) with three class labels (1, 2, 3) and only two features \(x_1\) and \(x_2\).
First, we define a latent variable as a linear combination of the features:
\[\begin{equation*} y_i^*=\beta_1 X_{i1}+\beta_2 X_{i2}+\mu_i \tag{6.63} \end{equation*}\]
Then since we have three classes we define two thresholds for this latent variable \(\alpha_1\) and \(\alpha_2\) (\(\alpha_1 < \alpha_2\)) such that a particular observation \(y_i\) will be classified as follows:
\[\begin{equation*} \begin{cases} y_i=1 & \text{if $y_i^* \leq \alpha_1$} \\ y_i=2 & \text{if $\alpha_1 < y_i^* \leq \alpha_2$} \\ y_i=3 & \text{if $y_i^* > \alpha_2$} \end{cases} \end{equation*}\]
Now we can obtain the probability of a particular observation falling into a specific class as follows (with non-stochastic regressors for simplicity to avoid conditionals):
\[\begin{equation} \begin{cases} p(y_i=1)=p(y_i^* \leq \alpha_1)=F(\alpha_1-\beta_1 X_{i1}-\beta_2 X_{i2}) \\ p(y_i=2)=p(\alpha_1 < y_i^* \leq \alpha_2)=F(\alpha_2-\beta_1 X_{i1}-\beta_2 X_{i2})-F(\alpha_1-\beta_1 X_{i1}-\beta_2 X_{i2}) \\ p(y_i=3)=1-p(y_i=2)-p(y_i=1) \end{cases} \tag{6.64} \end{equation}\]
Then, as discussed in 6.4, the suitable function can be the logistic function defined in 6.3.1 (\(\Gamma (x)=\frac{1}{1+exp^{-x}}\)). The above system then can be rewritten as follows:
\[\begin{equation*} \begin{cases} p(y_i=1)=\frac{1}{1+exp^{-(\alpha_1-\beta_1 X_{i1}-\beta_2 X_{i2})}} \\ p(y_i=2)=\frac{1}{1+exp^{-(\alpha_2-\beta_1 X_{i1}-\beta_2 X_{i2})}}-p(y_i=1) \\ p(y_i=3)=1-p(y_i=2)-p(y_i=1) \end{cases} \end{equation*}\]
We can further write that system in terms of the odds for interpretability. But let us treat each equation of the above system separately as follows (we set \(\gamma_1=\alpha_1-\beta_1 X_{i1}-\beta_2 X_{i2}\) and \(\gamma_2=\alpha_2-\beta_1 X_{i1}-\beta_2 X_{i2}\)):
The first equation:
\[\begin{align*} p(y_i=1)=\frac{1}{1+exp^{-\gamma_1}} & \implies p(y_i=1)+p(y_i=1)\times exp^{-\gamma_1} = 1 \\ & \implies \frac{p(y_i=1)}{1-p(y_i=1)}=exp^{-\gamma_1} \end{align*}\]
Since \(p(y_i=1)=p(y_i\leqslant 1)\) and \(1-p(y_i=1)=p(y_i>1)\), then by including the log, the equation will be:
\[\begin{equation*} log\Bigg[\frac{p(y_i\leqslant 1)}{p(y_i>1)}\Bigg]=\alpha_1-\beta_1 X_{i1}-\beta_2 X_{i2} \end{equation*}\]
The second equation:
Similarly, the second equation will be:
\[\begin{equation*} log\Bigg[\frac{p(y_i\leqslant 2)}{p(y_i>2)}\Bigg]=\alpha_2-\beta_1 X_{i1}-\beta_2 X_{i2} \end{equation*}\]
The probability of the last class can just be derived from the first once by complementarity as follows:
\[\begin{equation*} p(y_i=3)=1-p(y_i=2)-p(y_i=1) \end{equation*}\]
In general, if we have a set of regressors and ordinal response variable with \(k\) classes represented by integers \(j=1,..,k\), we will need \(k-1\) binary response models (since the last probability can be derived), each one will be defined by:
\[\begin{equation} log\Bigg[\frac{p(y\leqslant j)}{p(y>j)}\Bigg]=\alpha_j-X\beta \tag{6.65} \end{equation}\]
Any continuous distribution can be used instead of the logistic function. However, those implemented in Python and R are: probit, cloglog, loglog, cauchit, Aranda-Ordaz, log-gamma.
6.6.2 Continuation ratio model
It was proposed first by (Feinberg 1980). This model is similar to the last one with a slight difference which is instead of comparing the probability of being in one of the categories less or equal to the category (integer) \(j\) in the odds expression, we compare the probability of being beyond the category \(j\). Formally, this model is defined by:
\[\begin{equation} log\Bigg[\frac{p(y>j)}{p(y\geqslant j}\Bigg]=\alpha_j-X\beta \tag{6.66} \end{equation}\]
The regression coefficient \(\beta_i\) of a particular regressor \(x_i\) is the marginal effect on the log adds. Notice that the regression coefficients in \(\beta\) are identical across categories implying that every individual relationship between the response \(y\) and a particular \(x_i\) is independent of the classes.
6.6.3 Adjacent category logistic model
This model (Cande.V and David.G 1997) involves only the adjacent higher category. Formally, this model is defined by:
\[\begin{equation} log\Bigg[\frac{p(y=j)}{p(y=j+1)}\Bigg]=\alpha_j-X\beta_j \tag{6.67} \end{equation}\]
Notice that the regression coefficients \(\beta\)’s are allowed to change across categories in contrast to the above models.
Example 6.4 We will use the open source data diamonds available in kaggle:
Cumulative logit model:
import pandas as pd
# Read the data and remove the unnamed column
diam_py = pd.read_csv('diamonds.csv').iloc[:,1:]
diam_py.head(5)[out] cut color clarity depth table price x y z
[out] 0 Ideal E SI2 61.5 55.0 326 3.95 3.98 2.43
[out] 1 Premium E SI1 59.8 61.0 326 3.89 3.84 2.31
[out] 2 Good E VS1 56.9 65.0 327 4.05 4.07 2.31
[out] 3 Premium I VS2 62.4 58.0 334 4.20 4.23 2.63
[out] 4 Good J SI2 63.3 58.0 335 4.34 4.35 2.75The variable of our interest is cut, which is an ordinal variable with the following categories:
[out] cut
[out] Ideal 21551
[out] Premium 13791
[out] Very Good 12082
[out] Good 4906
[out] Fair 1610
[out] Name: count, dtype: int64As we see, this variable has five categories. To fit an ordinal model to this data in Python, we will call again the statsmodels package. However, we should first convert the type of that variable into a categorical one since it is now int64. This can be done very easily by using the pandas API module.
from pandas.api.types import CategoricalDtype
diam_py.cut=diam_py.cut.astype(CategoricalDtype(categories=["Fair","Good","Ideal","Very Good","Premium"], ordered=True))
# It's good practice to check what you expected
diam_py.cut.dtype[out] CategoricalDtype(categories=['Fair', 'Good', 'Ideal', 'Very Good', 'Premium'], ordered=True, categories_dtype=object)Now we are ready to fit an ordinal model using the cumulative logit model. However, for simplicity, we restrict the regressors set to include only the three last variables, the width x, the height y, and the depth z, then, we reduce the number of observations to \(10000\) instead of the actual number of \(53940\). We should also replace the small values in the variable \(z\) by the quantile \(0.01\)
diam_py = diam_py.loc[0:9999, ["cut","x", "y","z"]]
diam_py.loc[(diam_py.z < diam_py.z.quantile(0.01)), "z"] = diam_py.z.quantile(0.01)
diam_py.describe()[out] x y z
[out] count 10000.000000 10000.000000 10000.000000
[out] mean 5.985899 5.986773 3.703233
[out] std 0.652080 0.641055 0.406184
[out] min 3.790000 3.750000 2.500000
[out] 25% 5.780000 5.790000 3.560000
[out] 50% 6.130000 6.140000 3.810000
[out] 75% 6.400000 6.390000 3.970000
[out] max 7.620000 7.590000 4.870000To fit an ordinal model in Python, we call the function OrderedModel from statsmodels.
from statsmodels.miscmodels.ordinal_model import OrderedModel
mod_ordinal_py=OrderedModel(diam_py["cut"], diam_py[["x", "y","z"]], distr="logit")
res_ordinal_py=mod_ordinal_py.fit()
res_ordinal_py.params
[out] x 3.315428
[out] y 1.592731
[out] z -7.736140
[out] Fair/Good -2.480879
[out] Good/Ideal 0.495752
[out] Ideal/Very Good 0.454312
[out] Very Good/Premium 0.173076
[out] dtype: float64The first three parameters are the coefficients in the equation (6.63). However, the last three parameters are not equal exactly to the threshold values except for the first one (\(-2.480879\)). The others are exponentially incremented as follows:
import numpy as np
alpha_1=-2.480879
alpha_2=alpha_1+np.exp(0.495752)
alpha_3=alpha_2+np.exp(0.454312)
alpha_4=alpha_3+np.exp(0.173076)
print(f'alpha_1: {np.around(alpha_1, 4)}, alpha_2: {np.around(alpha_2, 4)}, alpha_3: {np.around(alpha_3, 4)}, alpha_4: {np.around(alpha_4, 4)} ')[out] alpha_1: -2.4809, alpha_2: -0.8391, alpha_3: 0.7359, alpha_4: 1.9249Or we can use the API as follows:
[out] array([ -inf, -2.48087853, -0.83914645, 0.73594295, 1.92489961,
[out] inf])If we plug in the coefficients and the computed thresholds in the above equations, we obtain the probability of each class as follows:
\[\begin{equation*} \begin{cases} p(y_i=Fair)=\frac{1}{1+exp^{-(-2.4808-3.3154x-1.5927y+7.7361z)}} \\ p(y_i=Good)=\frac{1}{1+exp^{-(-0.8391-3.3154x-1.5927y+7.7361z)}}-p(y_i=Fair) \\ p(y_i=Ideal)=\frac{1}{1+exp^{-(0.7359-3.3154x-1.5927y+7.7361z)}}-p(y_i=Fair)-p(y_i=Good) \\ p(y_i=Very Good)=\frac{1}{1+exp^{-(1.9249-3.3154x-1.5927y+7.7361z)}}-p(y_i=Fair)-p(y_i=Good)-p(y_i=Ideal) \\ p(y_i=Premium)=1-p(y_i=Very Good)-p(y_i=Fair)-p(y_i=Good)-p(y_i=Ideal) \end{cases} \end{equation*}\]
Now, if we predict the first observation, we should get the following probabilities:
x, y, z=diam_py[["x", "y","z"]].iloc[0]
fair=1/(1+np.exp(-(-2.4808-3.3154*x-1.5927*y+7.7361*z)))
good=(1/(1+np.exp(-(-0.8391-3.3154*x-1.5927*y+7.7361*z))))-fair
ideal=(1/(1+np.exp(-(0.7359-3.3154*x-1.5927*y+7.7361*z))))-fair-good
very_good=(1/(1+np.exp(-(1.9249-3.3154*x-1.5927*y+7.7361*z))))-fair-good-ideal
premium=1-fair-good-ideal-very_good
print(f'Fair:{np.around(fair,4)}, Good:{np.around(good,4)}, Ideal:{np.around(ideal,4)}, Very Good:{np.around(very_good, 4)}, Premium:{np.around(premium,4)}')[out] Fair:0.0707, Good:0.2114, Ideal:0.3729, Very Good:0.2067, Premium:0.1382However, we can get that prediction by merely calling the function predict.
# let us predict the first observation
predicted=mod_ordinal_py.predict(res_ordinal_py.params, exog=diam_py[["x", "y","z"]].values)
predicted[0][out] array([0.07072934, 0.21142356, 0.37289028, 0.20674696, 0.13820985])In R:
# convert the variable cut into an ordinal variable
diam_r$cut<-factor(diam_r$cut, levels=c("Fair", "Good", "Ideal", "Very Good", "Premium"), ordered = TRUE)
# using the first 10000 observations and the last three variables with the variable cut
diam_r <- diam_r[1:10000, c("cut","x","y","z")]
# correct the value of the variable z
diam_r$z <- ifelse(diam_r$z < quantile(diam_r$z, 0.01), quantile(diam_r$z, 0.01), diam_r$z )To fit a cumulative logistic model in R, we will use the R package ordinal.
[out] formula: cut ~ x + y + z
[out] data: diam_r
[out]
[out] link threshold nobs logLik AIC niter max.grad cond.H
[out] logit flexible 10000 -13943.84 27901.68 5(0) 3.97e-07 7.0e+04
[out]
[out] Coefficients:
[out] Estimate Std. Error z value Pr(>|z|)
[out] x 10.2314 0.3782 27.06 <2e-16 ***
[out] y -5.3003 0.3635 -14.58 <2e-16 ***
[out] z -7.9681 0.2315 -34.42 <2e-16 ***
[out] ---
[out] Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[out]
[out] Threshold coefficients:
[out] Estimate Std. Error z value
[out] Fair|Good -3.26425 0.17578 -18.570
[out] Good|Ideal -1.68629 0.16998 -9.920
[out] Ideal|Very Good -0.08115 0.16897 -0.480
[out] Very Good|Premium 1.15422 0.17017 6.783Using the above output, we can derive the probability of each class manually. Unlike Python, the threshold values are the same as in the output. For instance, to compute the probabilities of the first observation, we will do the following:
D=diam_r[1,c("x", "y","z")]
x <- D$x
y <- D$y
z <- D$z
fair=1/(1+exp(-(-3.2642-10.2314*x+5.3003*y+7.9681*z)))
good=1/(1+exp(-(-1.6863-10.2314*x+5.3003*y+7.9681*z)))-fair
ideal=1/(1+exp(-(-0.0811-10.2314*x+5.3003*y+7.9681*z)))-fair-good
very_good=1/(1+exp(-(1.1542-10.2314*x+5.3003*y+7.9681*z)))-fair-good-ideal
premium=1-fair-good-ideal-very_good
print(paste("Fair: ", round(fair, 4), "Good: ", round(good, 4), "Ideal: ", round(ideal, 4), "Very good: ", round(very_good, 4), "Premium: ", round(premium, 4)))
[out] [1] "Fair: 0.0652 Good: 0.1874 Ideal: 0.3746 Very good: 0.2254 Premium: 0.1473"However, we have computed those probabilities manually only for interpretation purposes because we can obtain them straightforwardly by calling the predict function as follows:
[out] $fit
[out] Fair Good Ideal Very Good Premium
[out] 1 0.06519378 0.1873642 0.3746201 0.2254662 0.1473558As we see the probabilities are the same.
Continuation ratio model:
In R:
To fit this model, we will use The vector generalized additive models VGAM package.
library(VGAM)
mod_ordinalc_r <- vglm(cut ~ x+y+z, cratio(parallel = TRUE), data = diam_r)
coef(mod_ordinalc_r, matrix = TRUE)[out] logitlink(P[Y>1|Y>=1]) logitlink(P[Y>2|Y>=2]) logitlink(P[Y>3|Y>=3])
[out] (Intercept) 3.075186 1.868450 0.348898
[out] x 10.255720 10.255720 10.255720
[out] y -6.136049 -6.136049 -6.136049
[out] z -6.628179 -6.628179 -6.628179
[out] logitlink(P[Y>4|Y>=4])
[out] (Intercept) -0.2469671
[out] x 10.2557199
[out] y -6.1360486
[out] z -6.6281785We have set parallel argument TRUE to obtain identical regression coefficients across the binary models. As we see, the coefficients are the same, but the thresholds are different
The prediction of the first observation using this model will be:
[out] P[Y>1|Y>=1] P[Y>2|Y>=2] P[Y>3|Y>=3] P[Y>4|Y>=4]
[out] 1 0.930433 0.8000527 0.4668168 0.3254596Since the actual value of the cut in the first observation is:
[out] [1] Ideal
[out] Levels: Fair < Good < Ideal < Very Good < PremiumThen, the two first probabilities are high because they match this level.
Adjacent category logistic model:
For this model, we will use the brglm2 package.
library(brglm2)
mod_ordinaladj_r <- bracl(cut ~ x+y+z, data = diam_r[1:1000,], parallel = TRUE, type = "ML")
summary(mod_ordinaladj_r)[out] Call:
[out] bracl(formula = cut ~ x + y + z, data = diam_r[1:1000, ], parallel = TRUE,
[out] type = "ML")
[out]
[out] Coefficients:
[out] Estimate Std. Error z value Pr(>|z|)
[out] Fair:(Intercept) -1.280 0.327 -3.92 8.9e-05 ***
[out] Good:(Intercept) -2.169 0.299 -7.25 4.1e-13 ***
[out] Ideal:(Intercept) -0.392 0.282 -1.39 0.1646
[out] Very Good:(Intercept) -0.948 0.280 -3.39 0.0007 ***
[out] x -3.721 0.582 -6.39 1.7e-10 ***
[out] y 2.204 0.545 4.05 5.2e-05 ***
[out] z 2.682 0.322 8.34 < 2e-16 ***
[out] ---
[out] Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[out]
[out] Residual Deviance: 2799.761
[out] Log-likelihood: -1399.88
[out] AIC: 2813.761
[out]
[out]
[out] Type of estimator: ML (maximum likelihood)
[out] Number of Fisher Scoring iterations: 4[out] Fair Good Ideal Very Good Premium
[out] 0.05272514 0.08677202 0.34737969 0.23523589 0.27788726In Python we will use the mord package as follows:
import pandas as pd
from mord import LogisticAT
diam_py['cut'] = diam_py['cut'].astype('category').cat.codes
X=diam_py[["x","y","z"]]
y = diam_py.cut
model = LogisticAT(alpha=1.0) # You can adjust the regularization parameter (alpha) as needed
model.fit(X, y)[out] LogisticAT()[out] array([[0.05944386, 0.15943151, 0.36284995, 0.24154479, 0.17672989]])