The classification problem of ordinal response data, examples of this including cancer patients grouped in early, mediocre and terminal stages, customers grouped into low, middle and high credit levels, are widely discussed in recent years. Some methods^[1-3] have been proposed to solve the problem. A direct method for ordinal response data is to convert ordinal labels to numerical values, such as the conversion of {excellent, good, moderate, poor} to {1, 2, 3, 4}. As a result, a regression method can be applied to the converted dataset. However, not all class labels can be converted to positive integers. Another method is to convert class labels to common numerical value rather than positive integers by a map. Nevertheless, it is difficult to define the mapping function. Frank and Hall^[4] proposed a method that converted ordinal classification problem to several binary classification problems. For each i∈{0, 1, 2, 3, …, K}, each sample is classified between the meta-class including classes 0 to i and the meta-class including classes i+1 to K. The final classification result can be inferred from the K binary predictions. Shashua and Levin^[5] extended SVM for ordinal regression by finding K ordinal thresholds, namely b₁≤b₂≤…≤b_K, spliting the real line into K+1 ordinal parts. Wang et al.^[6] proposed nonparallel support vector for ordinal regression by constructing a hyperplane in each class. However, these methods perform poorly when the ordinal label is mislabeled.

More and more literatures studied mislabeled response data by means of binary logistic regression^[7]. Tian and Sun^[8] proposed a new fuzzy set method to detect suspectable mislabeled points, and then delete their labels and construct a semi-supervised model. Logistic regression based on different mislabel probabilities are also widely used to mislabeled unordered response data. Copas considered equal and constant mislabel probabilities^[9] and Komori et al.^[10] assumed mislabeling occurs only the 0-group. In addition, Hung et al.^[7] proposed a robust mislabel logistic regression based on γ-divergence with the property of strong robustness in dealing with binary mislabeled response data. Neverthelss, above studies is not suitable to mislabeled ordinal response data.

In this article, we apply γ-divergence to ordinal multiple classification logistic regression and construct a robust ordinal mislabel logistic regression model, which not only possess the strong robustness but also the mislabel probabilities need not to be modeled. Our method have better effectiveness on ordinal mislabeled response data. It is noted that although γ-divergence has been previously adopted for regression analysis, to the best of our knowledge, γ-divergence has never been applied to ordinal response data.

1 Methodology

Given a dataset {(x_i, y_0i), i=1, 2, …, n} in binary logistic regression, where x_i∈ $\mathbb{R}$^p, y_0i∈{0, 1}, let p=P(y₀=1|X =x), 1-p=P(y₀=0|X =x), the log-odds (or logit) of probability p of an event is defined as logit(p)=log[p/(1-p)]=β^Tx+b.

1.1 Ordinal logistic regression

Suppose there are K+1 ordered classes {0, 1, 2, …, K}, π_j(x)=P(y₀=j|X=x). Consider the event that an observation belongs to the meta-class {0, 1, …, j}, that is y₀∈{0, 1, …, j}, whose probability is $\sum\limits_{k=0}^{j} \pi_{k}(\boldsymbol{x})$. The ordered logistic regression model (also called cumulative odds logit model^[11]) assumes that the log-odds are

$ \begin{aligned} &\operatorname{logit}\left[\pi_{0}(\boldsymbol{x})+\pi_{1}(\boldsymbol{x})+\pi_{2}(\boldsymbol{x})+\cdots+\pi_{j}(\boldsymbol{x})\right] \\ &=\log \left[\frac{\pi_{0}(\boldsymbol{x})+\pi_{1}(\boldsymbol{x})+\pi_{2}(\boldsymbol{x})+\cdots+\pi_{j}(\boldsymbol{x})}{\pi_{j+1}(\boldsymbol{x})+\pi_{j+2}(\boldsymbol{x})+\cdots+\pi_{K}(\boldsymbol{x})}\right] \\ &=\boldsymbol{\beta}^{\mathrm{T}} \boldsymbol{x}+b_{j+1}, j=0, \cdots, K-1 . \end{aligned} $

(1)

And the coefficient vector β=(β₁, β₂, …, β_p)^T is the same for all j in equation (1). Their intercept terms b₁, b₂, …, b_K split the plane into K+1 ordinal parts. Therefore, we only need to estimate parameters β, b₁, b₂, …, b_K, and β_j+1=(β^T, b_j+1)^T denotes (j+1)-th hyperplane, j=0, 1, 2, …, K-1.

According to Eq.(1), we obtain the category probability functions as follows:

$ \left\{\begin{array}{l} \pi_{0}\left(\boldsymbol{x} ; \boldsymbol{\beta}, b_{1}, \cdots, b_{K}\right)=\frac{A \exp \left(b_{1}\right)}{D} \\ \ \ \ \ \ \ \ \ \Pi_{h=2}^{K}\left[1+A \exp \left(b_{h}\right)\right], \\ \pi_{j}\left(\boldsymbol{x} ; \boldsymbol{\beta}, b_{1}, \cdots, b_{K}\right)=\frac{1}{D}\left[A \exp \left(b_{j+1}\right)-A \exp \left(b_{j}\right)\right] \\ \ \ \ \ \ \ \ \ \Pi_{l=1, l \neq j, l \neq j+1}^{K}\left[1+A \exp \left(b_{l}\right)\right], \\ \ \ \ \ \ \ \ \ j=1, \cdots, K-1, \\ \pi_{K}\left(\boldsymbol{x} ; \boldsymbol{\beta}, b_{1}, \cdots, b_{K}\right)=\frac{1}{D} \Pi_{h=1}^{K-1}\left[1+A \exp \left(b_{h}\right)\right], \end{array}\right. $

(2)

where $D=\Pi_{h=1}^{K}\left[1+A \exp \left(b_{h}\right)\right], A=\exp \left(\mathit{\pmb{\beta}}^{\mathrm{T}} \boldsymbol{x}\right)$.

1.2 The minimum γ-divergence estimation and its estimation equation

Let g be the data generating probability density function and f_θ be the model probability density function indexed by the parameter θ, the γ-divergence between g and f_θ is defined as^[7]

$ D_{\gamma}\left(g, f_{\boldsymbol{\theta}}\right)=\frac{1}{\gamma(\gamma+1)}\left\{\|g\|_{\gamma+1}-\int\left(\frac{f_{\boldsymbol{\theta}}}{\left\|f_{\boldsymbol{\theta}}\right\|_{\gamma+1}}\right)^{\gamma} g\right\}, $

(3)

where

$ \|g\|_{\gamma+1}=\left(\int g^{\gamma+1}\right)^{\frac{1}{\gamma+1}},\left\|f_{\boldsymbol{\theta}}\right\|_{\gamma+1}=\left(\int f_{\boldsymbol{\theta}}^{\gamma+1}\right)^{\frac{1}{\gamma+1}}. $

In the limiting case, $\lim\limits_{\gamma \rightarrow 0} D_{\gamma}\left(g, f_{\boldsymbol{\theta}}\right)=\int \ln \left(\frac{g}{f_{\boldsymbol{\theta}}}\right) g$, which is the KL-divergence^[7].

In the presence of contamination, g=εf_θ^*+(1-ε) τ, θ^* is the true model parameter, 1-ε is the contamination proportion and τ is the contamination density function. The estimation criterion of the minimum γ-divergence estimates parameter by minimizing D_γ(g, f_θ), which is equivalent to minimizing

$ \varepsilon D_{\gamma}\left(f_{\boldsymbol{\theta}^{*}}, f_{\boldsymbol{\theta}}\right)-\frac{H_{\gamma}(\varepsilon, \tau ; \boldsymbol{\theta})}{\gamma(\gamma+1)}, $

by taking the actual g into Eq.(3) and ignoring the terms which not involve θ, where H_γ(ε, τ; θ)=(1-ε) $\int\left(\frac{f_{\boldsymbol{\theta}}}{\left\|f_{\boldsymbol{\theta}}\right\|_{\gamma+1}}\right)^{\gamma} \boldsymbol{\tau}$.

Remark 1.1   Suppose that $\int f_{\theta^{*}}^{\gamma} \tau$ is sufficiently small for an appropriately large γ>0, which implies that the contamination density function τ mostly lies on the tail of the underlying density f_θ^*. Then for some γ, the bias H_γ(ε, τ; θ) is negligibly small when θ in a neighborhood of θ^*, this means that the estimation of θ is less effected by the 1-ε and τ. There are more introductions about γ-divergence and its property^[12] and more detailed discussions about the bias H_γ(ε, τ; θ)^[7].

1.3 Robust ordinal mislabel logistic regression method

Consider n independent samples and combine the discussions of H_γ(ε, τ; θ), the empirical version of the γ-divergence loss function^[7] is

$ L(\boldsymbol{\theta})=-\frac{1}{n} \sum\limits_{i=1}^{n} \frac{f\left(y_{i} \mid \boldsymbol{x}_{i} ; \boldsymbol{\theta}\right)^{\gamma}}{\left\|f\left(\cdot \mid \boldsymbol{x}_{i} ; \boldsymbol{\theta}\right)\right\|_{\gamma+1}}, $

(4)

where $\left\|f\left(\cdot \mid \boldsymbol{x}_{i} ; \boldsymbol{\theta}\right)\right\|_{\gamma+1}=\left(\int f\left(\mid \boldsymbol{x}_{i} ; \boldsymbol{\theta}\right)^{\gamma+1} \mathrm{~d} y\right)^{\frac{\gamma}{\gamma+1}}$.

In conventional ordinal logistic regression model with K+1 classes, namely y₀∈{0, 1, …, K}, the probability density function

$ f\left(y_{0} \mid \boldsymbol{x} ; \boldsymbol{\beta}, b_{1}, \cdots, b_{K}\right)=\prod\limits_{j=0}^{K} \pi_{j}^{Y_{0 j}}\left(\boldsymbol{x} ; \boldsymbol{\beta}, b_{1}, \cdots, b_{K}\right), $

(5)

for y₀, where π_j(x; β, b₁, …, b_K), j=0, 1, …, K given by the (2), and y₀ can be recorded by K random variables (Y₀₁, Y₀₂, …, Y_0K), where (0, 0, 0, …, 0) denotes y₀=0, y₀=j when Y_0j=1, Y_0k=0, j≠k. On the basis of (4), the robust estimator $\mathit{\pmb{\hat{{\beta}}}}_{\gamma}, \hat{b}_{j \gamma}, j=1, 2, \cdots, K$ can be acquired by minimizing the objective function

$ \begin{gathered} F\left(\boldsymbol{\beta}, b_{1}, \cdots, b_{K}\right)\\ =-\frac{1}{n} \sum\limits_{i=1}^{n}\left(\frac{f\left(y_{i} \mid \boldsymbol{x}_{i} ; \boldsymbol{\beta}, b_{1}, \cdots, b_{K}\right)}{\left\|f\left(\cdot \mid \boldsymbol{x}_{i} ; \boldsymbol{\beta}, b_{1}, \cdots, b_{K}\right)\right\|_{1+\gamma}}\right)^{\gamma}, \end{gathered} $

(6)

where

$ \begin{gathered} \left\|f\left(\cdot \mid \boldsymbol{x} ; \boldsymbol{\beta}, b_{1}, \cdots, b_{K}\right)\right\|_{1+\gamma} \\ =\left\{\sum\limits_{j=0}^{K}\left[\pi_{j}\left(\boldsymbol{x} ; \boldsymbol{\beta}, b_{1}, \cdots, b_{K}\right)\right]^{1+\gamma}\right\}^{\frac{1}{1+\gamma}} . \end{gathered} $

Direct differentiation of the objective function F(β, b₁, …, b_K) leads to the parameters estimation equations $S_{\gamma}\left(\mathit{\pmb{\hat{{\beta}}}}_{\gamma}\right)=0, S_{\gamma}\left(\hat{b}_{j \gamma}\right)=0, j=1, 2, \cdots, K$, where

$ S_{\gamma}(\boldsymbol{\beta})=-\frac{\gamma}{1+\gamma} \cdot \frac{1}{n} \sum\limits_{i=1}^{n} \omega_{\gamma i}\left(\frac{\frac{(1+\gamma) \partial C_{i}}{\partial \boldsymbol{\beta}}}{C_{i}}-\frac{\frac{\partial B_{i}}{\partial \boldsymbol{\beta}}} {B_{i}}\right), $

(7)

$ S_{\gamma}\left(b_{j}\right)=-\frac{\gamma}{1+\gamma} \cdot \frac{1}{n} \sum\limits_{i=1}^{n} \omega_{\gamma i}\left(\frac{\frac{(1+\gamma) \partial C_{i}}{\partial b_{j}}}{C_{i}}-\frac{\frac{\partial B_{i}}{\partial b_{j}}}{B_{i}}\right), $

(8)

with the weight function

$ {\omega}_{\gamma i}=\left(\frac{C_{i}^{1+\gamma}}{{B}_{i}}\right)^{\frac{\gamma}{\gamma+1}}, $

(9)

and where $A_{i}=\exp \left(\boldsymbol{\beta}^{\mathrm{T}} \boldsymbol{x}_{i}\right)$,

$ \begin{aligned} &B_{i}=\\ &\sum\limits_{j=1}^{K-1}\left\{\left[A_{i}\left(\exp \left(b_{j+1}\right)-\exp \left(b_{j}\right)\right)\right] \prod\limits_{l=1,\atop {l \neq j,\atop l \neq j+1}}^{K}\left[1+A_{i} \exp \left(b_{l}\right)\right]\right\}^{1+\gamma}+\\ &\left[A_{i} \exp \left(b_{1}\right) \prod\limits_{h=2}^{K}\left(1+A_{i} \exp \left(b_{h}\right)\right)\right]^{1+\gamma}+\prod\limits_{j=1}^{K-1}\left[1+A_{i} \exp \left(b_{j}\right)\right]^{1+\gamma}, \end{aligned} $

(10)

$ \begin{aligned} &C_{i}= \\ &\sum\limits_{j=1}^{K-1}\left\{Y_{i j}\left[A_{i}\left(\exp \left(b_{j+1}\right)-\exp \left(b_{j}\right)\right)\right] \prod\limits_{\substack{l=1, \atop {l \neq j,\atop l \neq j+1}}}^{K}\left[1+A_{i} \exp \left(b_{l}\right)\right]\right\}+ \\ &\left(1-\sum\limits_{j=1}^{K} Y_{i j}\right) A_{i} \exp \left(b_{1}\right) \prod\limits_{h=2}^{K-1}\left[1+A_{i} \exp \left(b_{h}\right)\right]+ \\ &Y_{i K} \prod\limits_{h=1}^{K-1}\left[1+A_{i} \exp \left(b_{h}\right)\right] \end{aligned} $

(11)

When γ=0, the above estimation equation degenerates to the estimation equation of conventional ordinal logistic regression with non-robust estimator. From (7)-(9), the robustness of $\mathit{\pmb{\hat{{\beta}}}}_{\gamma}, \hat{b}_{j \gamma}, j=1, 2, \cdots, K$ is clear, that's to say, when mislabeled, the weight ω_γi of the sample x_i diminish in comparison with the case that the x_i is properly labeled. This point is explained at the Remark 1.2.

Remark 1.2   The weight of x in the parameter estimation equation is $\omega_{\gamma}=\left(\frac{C^{1+\gamma}}{B}\right)^{\frac{\gamma}{1+\gamma}}$. Let the ω_{γ, k} represents the weight of x when it is judged as the k-th class. When x is mislabeled, the weight diminish in comparison with the situation that the sample x is properly labeled, which can be described as the following proposition: if the true label of the sample x is j, i.e. π_j>π_h, then ω_{γ, j}>ω_{γ, h}, j=0, 1, 2, …, K, h≠j. For example, when K=2, we assume that the true label of the x is 1, we demonstrate π₁>π₀, π₁>π₂$\Rightarrow$ω_{γ, 1}>ω_{γ, 0}, ω_{γ, 1}>ω_{γ, 2} correspondingly. $\pi_{0}=\frac{A \exp \left(b_{1}\right)}{1+A \exp \left(b_{1}\right)}$, π₁= $\frac{A \exp \left(b_{2}\right)}{1+A \exp \left(b_{2}\right)}-\frac{A \exp \left(b_{1}\right)}{1+A \exp \left(b_{1}\right)}$, we get exp(b₂)>2exp(b₁)+Aexp(b₁+b₂) when π₁>π₀. ω_{γ, 0}= $\left\{\frac{\left[A \exp \left(b_{1}\right)\left(1+A \exp \left(b_{2}\right)\right)\right]^{1+\gamma}}{B}\right\}^{\frac{\gamma}{1+\gamma}}$, ω_{γ, 1}= $\left\{\frac{\left[A \exp \left(b_{2}\right)-A \exp \left(b_{1}\right)\right]^{1+\gamma}}{B}\right\}^{\frac{\gamma}{1+\gamma}}$, if ω_{γ, 1}>ω_{γ, 0} we need exp(b₂)>2exp(b₁)+Aexp(b₁+b₂), so π₁>π₀$\begin{equation}\Rightarrow\end{equation}$ω_{γ, 1}>ω_{γ, 0}. Similarly, when π₁>π₂, we can deduce ω_{γ, 1}>ω_{γ, 2}. Similar conclusions can be drawn under other K.

1.4 Computation

Due to the complexity of the second derivative of the objective function based on γ-divergence, in this paper, we adopt the gradient descent algorithm to solve model, which is summarized in Algorithm 1.

Algorithm 1   Gradient descent algorithm
Require: The dataset: (xi, yi), i=1, 2, …, n, the tuning parameter: γ.
Ensure: The regression coefficients: β, b1, b2, …, bK. Initializing the s=0, ε=10-3 and given an initial value β(0), bj(0), j=1, 2, …, K.
repeat
s=s+1;
for k=1, 2, …, p do
$\frac{\partial F}{\partial \beta_{k}}=-\frac{\gamma}{1+\gamma} \cdot \frac{1}{n} \sum\limits_{i=1}^{n} \omega_{\gamma i}\left( {\frac{{\frac{{(1 + \gamma )\partial {C_i}}}{{\partial {\beta _k}}}}}{{{C_i}}} - \frac{{\frac{{\partial {B_i}}}{{\partial {\beta _k}}}}}{{{B_i}}}} \right);$
end for
for j=1, 2, …, K do
$\frac{\partial F}{\partial b_{j}}=-\frac{\gamma}{1+\gamma} \cdot \frac{1}{n} \sum\limits_{i=1}^{n} \omega_{\gamma i}\left( {\frac{{\frac{{(1 + \gamma )\partial {C_i}}}{{\partial {b_j}}}}}{{{C_i}}} - \frac{{\frac{{\partial {B_i}}}{{\partial {b_j}}}}}{{{B_i}}}} \right);$
end for
$\begin{gathered} \left(\boldsymbol{\beta}^{(s+1)}, \boldsymbol{b}^{(s+1)}\right)=\left(\boldsymbol{\beta}^{(s)}, \boldsymbol{b}^{(s)}\right)- \\ \left.\alpha \frac{\partial F}{\partial(\boldsymbol{\beta}, \boldsymbol{b})}\right|_{\boldsymbol{\beta}=\boldsymbol{\beta}^{(s)}, b=b^{(s);}} \end{gathered}$
Where β=(β1, β2, …, βp)T, b=(b1, b2, …, bK)T, ωγi, Bi, Ci are defined as (9)-(11), α is the step length obtained by Armijo search. We consider three classification problem in our simulation study, the calculation results of $\begin{equation} \frac{\partial C_{i}}{\partial \beta_{k}}, \frac{\partial B_{i}}{\partial \beta_{k}}, \frac{\partial C_{i}}{\partial b_{j}}, \frac{\partial B_{i}}{\partial b_{j}} \end{equation}$ are placed in Appendix (A1)-(A6).
until
$\frac{\left.\|\left(\boldsymbol{\beta}^{\mathrm{T}(\mathrm{s}+1)}\right), \boldsymbol{b}^{\mathrm{T}(\mathrm{s}+1)}\right)-\left(\boldsymbol{\beta}^{\mathrm{T}(\mathrm{s})}, \boldsymbol{b}^{\mathrm{T}(\mathrm{s})}\right) \|_{2}}{\|\left(\boldsymbol{\beta}^{\mathrm{T}(\mathrm{s})}, \boldsymbol{b}^{\mathrm{T}(\mathrm{s})} \|_{2}\right.} \leqslant \boldsymbol{\varepsilon} .$
returnβT(s+1), bT(s+1).

The proposed objective function has a complex form, so it is difficult to establish the convergence property. Experimental results demonstrate that, for all of our simulated and real datasets, convergence is successfully achieved within 50 overall iterations (mostly within 20 iterations).

The definition of γ-divergence and studies suggest that γ balances between robustness and efficiency. However, Fujisawa and Eguchi^[12] said there could be not consistent best way to select an appropriate tuning parameter γ. In practice, there are some methods to select the tuning parameter γ, such as the adaptive selection procedure by Mollah et al.^[13], the cross validation^[14] and a sequence of the parameter^[15]. In this article, we consider a sequence of γ as 1, 2, 3, 4 and 5 following Zang et al.^[15], and proceed the simulation study under each determined γ.

2 Simulation studies

In the numerical simulation study, we consider ordinal three classification problem with possibly mislabeling, that's to say, the labels of ordinal response datas y∈{0, 1, 2}, then π₀(x), π₁(x), π₂(x), the probability density function f(y|x; β, b₁, b₂), and C_i, B_i, A_i can be obtained by the (2), (5), (10)-(11), and the calculation results of them are placed in Appendix (A7)-(A13). Simulation results are reported with 500 replicates in the simulation studies.

2.1 Simulation settings

In each simulation run, we generate 300 random samples x_s in $\mathbb{R} $¹⁴, and x_s~N(0, Σ), we consider three different structures of covariance matrix Σ=(σ_ij)_{1≤i, j≤14},

1) the explanatory variables are independent, i.e. Σ=I,

2) the auto-regressive correlation (AR) given by σ_ij=ρ^{|i-j |} with ρ=0.5 and ρ=0.75,

3) the banded correlation (Band) given by σ_ij=ρ^{|i-j |} with ρ=0.6 if |i-j|≤2, ρ=0 if |i-j|>2.

We generate ordered mislabeled labels {0, 1, 2} in the light of the following settings of mislabel probabilities:

$ \left\{\begin{array}{l} \eta_{01}=\eta_{10}=m / 2, \\ \eta_{12}=\eta_{21}=m / 2, \\ \eta_{02}=\eta_{20}=m / 4, \end{array}\right. $

(12)

where

$ \eta_{i j}=P\left(y=j \mid y_{0}=i, \boldsymbol{X}=\boldsymbol{x}\right), i, j=0,1,2. $

We assume m={0.1, 0.15, 0.2, 0.25, 0.3} in the simulation. Considering the mislabel probabilities, the labels can be generated by following formulas:

$ \left\{\begin{array}{l} P(y=0 \mid \boldsymbol{X}=\boldsymbol{x})=\left(1-\eta_{01}-\eta_{02}\right) \pi_{0}+\eta_{10} \pi_{1}+\eta_{20} \pi_{2}, \\ P(y=1 \mid \boldsymbol{X}=\boldsymbol{x})=\left(1-\eta_{10}-\eta_{12}\right) \pi_{1}+\eta_{01} \pi_{0}+\eta_{21} \pi_{2}, \\ P(y=2 \mid \boldsymbol{X}=\boldsymbol{x})=\left(1-\eta_{20}-\eta_{21}\right) \pi_{2}+\eta_{02} \pi_{0}+\eta_{12} \pi_{1} . \end{array}\right. $

(13)

2.2 Simulation results

We apply the robust ordinal mislabel logistic regression method on the simulated dataset. And we compare the conventional ordinal logistic regression (COLR) method which takes no account of mislabeling with the robust ordinal mislabel logistic regression (robust ordinal mislabel) method from different aspects. Firstly, we evaluate the two methods via two indexes, i.e. the mean absolute error (MAE) and standard deviation (SD) of the estimates of β, b₁, b₂. Secondly, we compare the classification accuracy of the two methods on 800 test samples. The true values, namely β₀, b₀₁, b₀₂ follow the first part of uniform distribution [-1, 0.5] and the last part of uniform distribution [0.5, 1], and β₀₁=(β₀^T, b₀₁)^T, β₀₂=(β₀^T, b₀₂)^T are two parallel classification hyperplanes. Partial results are presented in the article.

Figure 1 displays the variation of classification error rate of 800 test samples with m under different structures of Σ. It can be seen clearly that the classification error rate obtained by COLR method increases more greatly than the classification error rate acquired by robust ordinal mislabel method with the increase of the mislabel probability (MP). Furthermore, we can select γ by the performance of classification from Fig. 1. We can also note that there exists the best effectiveness when γ=4.

When Σ=I, the comparisons between the true values and the estimates of β, b₁, b₂ under different γ and m are presented in Fig. 2 (a)-Fig. 2 (d) and the comparisons between the true values and the estimates of β, b₁, b₂ under different γ and m are displayed in Fig. 2 (e)-Fig. 2 (h) when the structure of Σ is AR with ρ=0.5. Figure 2 show that, for all given m and γ and the structure of Σ is AR (ρ=0.5) or I, the estimates of β, b₁, b₂ from the robust ordinal mislabel method are closer to the true values of β, b₁, b₂ than from COLR method.

When the structure of Σ is Band, MAE of parameter estimation under different γ and m values are displayed in Fig. 3 (a)-Fig. 3 (d). When the structure of Σ is AR with ρ=0.75, the MAE of parameter estimation are presented in Fig. 3 (e)-Fig. 3 (h). Figure 3 illustrate that the MAE of parameter estimation from the robust ordinal mislabel method is smaller than from COLR method. Both Fig. 2 and Fig. 3 demonstrate the robust ordinal mislabel method possesses better performance in estimating coefficients of the classification hyperplanes.

Table 1 displays MAE of parameter estimation and SD of the estimates from the COLR method and the robust ordinal mislabel method under different Σ, MP and γ values. We can obviously find that the MAE and the SD from robust ordinal mislabel method are smaller than from COLR method respectively. It also demonstrates that the robust ordinal mislabel method performs well and more robustly than COLR method in estimating coefficients.

In addition, we verify the effectiveness of the proposed method on samples conforming to the heavy tailed distribution, such as t-distribution, mixture of normal distribution. In each simulation, we generate two groups of random samples x_s in $\mathbb{R}$¹⁴, 300 in each group. x_s are from the distribution t(0, Σ, df=3) in the first group and x_s are from mixture of normal distribution (0.2, 0, I; 0.3, 0, AR(ρ=0.5); 0.3, 0, AR(ρ=0.75); 0.2, 0, Band) in the second group.

When γ=4 and x_s are from the mixture of normal distribution, the MAE and SD of parameter estimation are displayed in Table 2. When γ=2 and x_s are from the distribution t(0, Σ, 3), the MAE and SD of parameter estimation are displayed in Table 3. From the resluts, we can obviously find that both MAE and SD from the COLR mtehod are larger than both MAE and SD from the robust ordinal mislabel method when x_s are from heavy tailed distributions. This also means that the robust ordinal mislabel method has better effectiveness than COLR method in estimating coefficients. Simultaneously, we can find from the Table 2 and Table 3 that MAE of parameter estimation is smaller when samples are from normal distribution in contrast with the situation that samples are from t-distribution or mixture of normal distribution.

3 Real data analysis

In this section, we demonstrate the effectiveness of the proposed method on real data. Our dataset, i.e. the childhood eye data, is from Beijng Tongren Hospital^[16], which contains 6 years of data from the first grade to sixth grade, and the variables include ocular biometry, near work, food habits, living habits, habits of wearing spectacles in this school year, accommodative response, times outdoors, and parental myopia and so on.

In the childhood vision research, we regard the change of spherical equivalent (CSE) after mydriasis from first grade to sixth grade as dependent variable. We divide the range of myopia into high myopia, moderate myopia and low myopia. The CSE < -2.5 denotes high myopia, corresponding to class 2 and the label is 2 in this class. The -2.5≤ CSE≤-1.5 denotes moderate myopia and corresponds to class 1. The CSE>-1.5 denotes low myopia and corresponds to class 0. The level of myopia is an ordinal response data and spherical equivalent might be inaccurate because of limitation of mydriasis examination, so there might be some mislabeled responses in the data.

We use the factors of vision from the first grade students as the independent variables to construct regression model, which can analyze the effectivenss of variables on the myopia progression in primary school. These variables are divided into three categories: continuous variables, nominal variables and multi-class variables. The data used in the model contains 1 370 samples and 28 independent variables and the meanings of the variables are presented in Appendix B.

We analyse the real data from two cases.

Example 3.1   (assessment of robustness) We first delete 10 % samples before running and then calculate the mean (Mean) and SD of ten runs of each variable. The Mean and SD are presented in Table 4. Our method performs more robustly than COLR by SD obviously. Most important variables^[17], such as "AL", "SE", "DUCVA", "D_COMR2", "BREAK3", have smaller SD by means of robust ordinal mislabel method than COLR method. It demonstrates that our method have strong stability.

Example 3.2  We take the Mean as the estimates of β, b₁, b₂. The absolute value of the Mean can reflect the importance of variables. Our method can better reflect the importance of important variables, such as "AL", "SE", which are regarded as important variables^[18], having greater absolute estimate with the help of robust ordinal mislabel method than COLR method. "LT" is positively correlated with the degree of myopia^[19]. However, the estimate of it from robust ordinal mislabel method is positive than from COLR method is negative. It shows that the robust ordinal mislabel method proposed in the paper can more accurately reflect the relationship between independent variables and dependent variables. In addition, "K1" is negatively correlated with the degree of myopia. It can reflect this point from the case of γ=1, 4, 5. However, the case of γ=2, 3 can not reflect the situation. The absolute value of the estimate of "K1" is the largest for the case of γ=4 and the result is consistent with classification error rate under different γ.

4 Conclusion and future work

In this paper, we proposed a robust ordinal mislabel logistic regression method based on γ-divergence. The model is obtained by minimizing γ-divergence estimation. Both theoretical analysis and simulation studies make clear that the proposed method is quite efficient for solving classification problems with ordinal and mislabeled response data. Firstly, the mislabel probabilities need not to be modeled in our method. Secondly, our method performs more robustly than conventional ordinal logistic regression by the ways of simulation results and real data analysis. However, the robust ordinal mislabel logistic regression based on γ-divergence is applied to low dimensional data. There exists mislabeled response data on some high dimensional ordinal data, we can consider further research for ordinal mislabel method to high dimensional data.

Appendices Appendix A

We consider ordinal three classification problem in simulation studies, the $\begin{equation} \frac{\partial C_{i}}{\partial \beta_{k}}, \frac{\partial B_{i}}{\partial \beta_{k}}, \frac{\partial C_{i}}{\partial b_{j}}, \frac{\partial B_{i}}{\partial b_{j}} \end{equation}$ in Algorithm 1 are as follows,

$ \frac{\partial C_{i}}{\partial b_{1}}=\left(1-Y_{i 1}-Y_{i 2}\right) A_{i} \exp \left(b_{1}\right)\left(1+A_{i} \exp \left(b_{2}\right)\right)-Y_{i 1} A_{i} \exp \left(b_{1}\right)+Y_{i 2} A_{i} \exp \left(b_{2}\right). $

(A1)

$ \begin{gathered} \frac{\partial B_{i}}{\partial b_{1}}=\left\{-\left[A_{i}\left(\exp \left(b_{2}\right)-\exp \left(b_{1}\right)\right)\right]^{\gamma} A_{i} \exp \left(b_{1}\right)+\left(1+A_{i} \exp \left(b_{1}\right)\right)^{\gamma} A_{i} \exp \left(b_{1}\right)+\right. \\ \left.\left[A_{i} \exp \left(b_{1}\right)\left(1+A_{i} \exp \left(b_{2}\right)\right)\right]^{1+\gamma}\right\}(1+\gamma) . \end{gathered} $

(A2)

$ \frac{\partial C_{i}}{\partial b_{2}}=\left(1-Y_{i 1}-Y_{i 2}\right) A_{i}^{2} \exp \left(b_{1}+b_{2}\right)+Y_{i 1} A_{i} \exp \left(b_{2}\right) . $

(A3)

$ \frac{\partial B_{i}}{\partial b_{2}}=\left\{\left[A_{i}\left(\exp \left(b_{2}\right)-\exp \left(b_{1}\right)\right)\right]^{\gamma} A_{i} \exp \left(b_{2}\right)+\left[A_{i} \exp \left(b_{1}\right)\left(1+A_{i} \exp \left(b_{2}\right)\right)\right]^{\gamma} A_{i}^{2} \exp \left(b_{1}+b_{2}\right)\right\}(\gamma+1). $

(A4)

$ \begin{gathered} \frac{\partial C_{i}}{\partial \beta_{k}}=\left\{Y_{i 1}\left[A_{i}\left(\exp \left(b_{2}\right)-\exp \left(b_{1}\right)\right)\right]+Y_{i 2} A_{i} \exp \left(b_{1}\right)+\right. \\ \left.\left(1-Y_{i 1}-Y_{i 2}\right)\left[A_{i}^{2} \exp \left(b_{1}+b_{2}\right)+A_{i} \exp \left(b_{1}\right)\left(1+A_{i} \exp \left(b_{2}\right)\right)\right]\right\} x_{i j} \cdot \end{gathered} $

(A5)

$ \begin{gathered} \frac{\partial B_{i}}{\partial \beta_{k}}=(1+\gamma)\left\{\left[A_{i}\left(\exp \left(b_{2}\right)-\exp \left(b_{1}\right)\right)\right]^{\gamma+1}+A_{i} \exp \left(b_{1}\right)\left(1+A_{i} \exp \left(b_{1}\right)\right)^{\gamma}+\right. \\ \left.\left[A_{i} \exp \left(b_{1}\right)\left(1+A_{i} \exp \left(b_{2}\right)\right)\right]^{1+\gamma}+\left[A_{i} \exp \left(b_{1}\right)\left(1+A_{i} \exp \left(b_{2}\right)\right)\right]^{\gamma} A_{i} \exp \left(b_{1}+b_{2}\right)\right\} x_{i j}. \end{gathered} $

(A6)

$ \pi_{0}=P\left(y_{0}=0 \mid \boldsymbol{X}=\boldsymbol{x}\right)=\frac{A \exp \left(b_{1}\right)}{1+A \exp \left(b_{1}\right)}. $

(A7)

$ \pi_{1}=P\left(y_{0}=1 \mid \boldsymbol{X}=\boldsymbol{x}\right)=\frac{A \exp \left(b_{2}\right)}{1+A \exp \left(b_{2}\right)}-\frac{A \exp \left(b_{1}\right)}{1+A \exp \left(b_{1}\right)} . $

(A8)

$ \pi_{2}=P\left(y_{2}=2 \mid \boldsymbol{X}=\boldsymbol{x}\right)=\frac{1}{1+A \exp \left(b_{2}\right)}. $

(A9)

The objective function:

$ F\left(\boldsymbol{\beta}, b_{1}, b_{2}\right)=-\frac{1}{n} \sum\limits_{i=1}^{n}\left(\frac{f\left(y_{i} \mid \boldsymbol{x}_{i} ; \boldsymbol{\beta}, b_{1}, b_{2}\right)}{\left\|f\left(\cdot \mid \boldsymbol{x}_{i} ; \boldsymbol{\beta}, b_{1}, b_{2}\right)\right\|}_{1+\gamma}\right)^{\gamma}=-\frac{1}{n} \sum\limits_{i=1}^{n}\left(\frac{C_{i}^{1+\gamma}}{B_{i}}\right)^{\frac{\gamma}{1+\gamma}}, $

(A10)

where $\begin{equation} \left\|f\left(\cdot \mid \boldsymbol{x} ; \mathit{\pmb{β}}, b_{1}, b_{2}\right)\right\|_{1+\gamma}=\left\{\sum_{j=0}^{2} \pi_{j}^{1+\gamma}\right\}^{\frac{1}{1+\gamma}} \end{equation}$,

$ A_{i}=\exp \left(\boldsymbol{\beta}^{\mathrm{T}} \boldsymbol{x}_{i}\right), $

(A11)

$ B_{i}=\left\{A_{i} \exp \left(b_{1}\right)\left[1+A_{i} \exp \left(b_{2}\right)\right]\right\}^{1+\gamma}+\left[A_{i}\left(\exp \left(b_{2}\right)-\exp \left(b_{1}\right)\right)\right]^{1+\gamma}+\left[1+A_{i} \exp \left(b_{1}\right)\right]^{1+\gamma}, $

(A12)

$ C_{i}=\left(1-Y_{i 1}-Y_{i 2}\right) A_{i} \exp \left(b_{1}\right)\left[1+A_{i} \exp \left(b_{2}\right)\right]+Y_{i 1}\left[A_{i}\left(\exp \left(b_{2}\right)-\exp \left(b_{1}\right)\right)\right]+Y_{i 2}\left[1+A_{i} \exp \left(b_{1}\right)\right]. $

(A13)

Appendix B

       The meanings of the variables in Table 4 are as follows. The level number after some variable name is the result of the dummy variable processing.

       Binary variables (1: yes, 0: no):

MOTHER_A: Whether the natural mother has amblyopia?

GLASS: Do you wear glasses?

DESK_L: When you are reading or doing close work, do you use a desk lamp?

OMT: Do your children have received the other myopia treatment?

TUTOR2: Do your children take part in the tutoring class in their spare time (indoor learning class)?

GENDER: Is gender a woman?

       Ordinal variables (1:lower, 4:upper):

COLA: Quantity of drinking a carbonated drink frequency in the last 4 weeks.

EGGEAT: Quantity of eating eggs frequency in the last 4 weeks.

MEATS: Quantity of eating meat frequency in the last 4 weeks.

READLY: Quantity of weekly reading.

BREAK: Quantity of keeping reading or doing close work frequency.

D_COMR: When your child uses a computer, distance from computer from the computer (1:nearer, 4:farer).

       Continuous variables:

MYOPICS: The number of myopia parents.	K2: The steep keratometry readings.
SE: Spherical equivalent after mydriasis.	AR: Accommodative response.
DUCVA: Distant uncorrected visual acuity.	PUPIL_SIZE: Pupil diameter.
AL: Axial length.	BL: Birth length.
LT: Lens thickness.	NUCVA: Near uncorrected visual acuity.
K1: The flat keratometry readings.