Parameter Estimation and Hypothesis Testing of The Bivariate Polynomial Ordinal Logistic Regression Model

Abstract

Logistic regression is one of statistical methods that used to analyze the correlation between categorical response variables and predictor variables which are categorical or continuous. Many studies on logistic regression have been carried out by assuming that the predictor variable and its logit link function have a linear relationship. However, in several cases it was found that the relationship was not always linear, but could be quadratic, cubic, or in the form of other curves, so that the assumption of linearity was incorrect. Therefore, this study will develop a bivariate polynomial ordinal logistic regression (BPOLR) model which is an extension of ordinal logistic regression, with two correlated response variables in which the relationship between the continuous predictor variable and its logit is modeled as a polynomial form. There are commonly two correlated response variables in scientific research. In this study, each response variable used consisted of three categories. This study aims to obtain parameter estimators of the BPOLR model using the maximum likelihood estimation (MLE) method, obtain test statistics of parameters using the maximum likelihood ratio test (MLRT) method, and obtain algorithms of estimating and hypothesis testing for parameters of the BPOLR model. The results of the first partial derivatives are not closed-form, thus, a numerical optimization such as the Berndt–Hall–Hall–Hausman (BHHH) method is needed to obtain the maximum likelihood estimator. The distribution statistically test is followed the Chi-square limit distribution, asymptotically.

1. Introduction

The logistic regression model is one of statistical methods which is used to analyze the correlation between categorical response variables and predictor variables that are categorical or continuous. If the response variable has more than two categories and there are levels in that category (ordinal scale), then an ordinal logistic regression model is used [1]. The logistic regression modeling often does not involve only one response variable because of a phenomenon involving multiple response variables. Logistic regression involving one response variable is called univariate logistic regression. Meanwhile, logistic regression that involves two or more response variables and correlates between response variables is called multivariate logistic regression. Especially for multivariate logistic regression involving two correlated response variables, it is called the bivariate logistic regression.

Several studies on ordinal logistic regression have been carried out for bivariate cases. Dale [2] has studied parameter estimationusing the MLE method. Williamson et al. [3] used the generalized estimating equation (GEE) method, whereas Enea and Lovison [4] used penalized MLE for the parameter estimation of bivariate ordinal logistic regression models. The previous studies examined logistic regression modeling by assuming that the continuous predictor variable and its logit link function have a linear relationship. However, the reality is that in some cases it is found that the relationship is not always linear, but can be quadratic, cubic, or other curves, so that the assumption of linearity is incorrect.

Royston and Altman [5] introduced the fractional polynomial approach, which is one of the development of a polynomial model, to understand the functional form of a continuous predictor variable. This model involves a more flexible relationship that can be linear or non-linear and has been used by several researchers [6,7,8,9]. In its development, several studies have examined the use of fractional polynomials in logistic regression analysis [10,11,12]. However, these previous studies were limited to univariate binary logistic regression. Furthermore, the bivariate binary logistic regression has been developed but limited to the second order [13]. Therefore, we propose an extension of ordinal logistic regression with two correlated response variables in which the relationship between the continuous predictor variable and its logit is modeled as a polynomial form. The proposed model is called the bivariate polynomial ordinal logistic regression (BPOLR) model. This study aimed to obtain parameter estimators of the BPOLR model using the MLE method, obtain test statistics of parameters using the MLRT method, and obtain algorithms of estimating and hypothesis testing for parameters of the BPOLR model.

An example of applying the BPOLR model can be applied to modeling the nutritional status of toddlers. Based on Narendra et al. [14], children from birth to about one year old have growth that increases rapidly and then decreases slowly as the child gets older. Therefore, the polynomial approach is more suitable for modeling child growth [15]. Tilling et al. [16] modeled child growth using the fractional polynomial approach and concluded that the model with the fractional polynomial approach shows a growth rate that is initially fast and then slows down over time.

The following discussion in this paper is divided into several main topics. In Section 2, we introduce the BPOLR model followed by parameter estimation using the MLE method in Section 3. We present simultaneous hypothesis testing in Section 4 to test the significance of all parameters together using the MLRT method and we present a simulation study in Section 5. Furthermore, Section 6 contains conclusions.

2. The Bivariate Polynomial Ordinal Logistic Regression (BPOLR) Model

The BPOLR model is an expansion of the ordinal logistic regression model when there are two correlated ordinal response variables and the relationship between the continuous predictor variable and its logit is modeled as a polynomial form. In this study, the response variables used each have 3 categories. Let Y₁ and Y₂ be response variables which have a value of 1; 2; or 3, then ${\mathrm{Y}}_{\mathrm{ab}}=\left({\mathrm{Y}}_1{=\mathrm{a},\mathrm{Y}}_2=\mathrm{b}\right) ; \mathrm{a},\mathrm{b}=1,2,3$ are random variables that have their respective probabilities ${\mathsf{\pi }}_{\mathrm{ab} }$, as presented in

Table 1. Probabilty of the response variables.

${\mathbf{Y}}_1$	${\mathbf{Y}}_2$			Total
	1	2	3
1	${\mathsf{\pi }}_{11}$	${\mathsf{\pi }}_{12}$	${\mathsf{\pi }}_{13}$	${\mathsf{\pi }}_{1•}$
2	${\mathsf{\pi }}_{21}$	${\mathsf{\pi }}_{22}$	${\mathsf{\pi }}_{23}$	${\mathsf{\pi }}_{2•}$
3	${\mathsf{\pi }}_{31}$	${\mathsf{\pi }}_{32}$	${\mathsf{\pi }}_{33}$	${1- \mathsf{\pi }}_{1•}{- \mathsf{\pi }}_{2•}$
Total	${\mathsf{\pi }}_{•1}$	${\mathsf{\pi }}_{•2}$	${1- \mathsf{\pi }}_{•1}{- \mathsf{\pi }}_{•2}$	1

below.

Based on Table 1, the random vector $\mathbf{y}=\left[\begin{array}{cccccccc}{\mathrm{Y}}_{11}& {\mathrm{Y}}_{12}& {\mathrm{Y}}_{13}& {\mathrm{Y}}_{21}& {\mathrm{Y}}_{22}& {\mathrm{Y}}_{23}& {\mathrm{Y}}_{31}& {\mathrm{Y}}_{32}\end{array}\right]$ have a multinomial distribution so it has a joint probability density function as follows:

$$\mathrm{P}\left({\mathrm{Y}}_{11}{=\mathrm{y}}_{11}{,\mathrm{Y}}_{12}{=\mathrm{y}}_{12},\dots {,\mathrm{Y}}_{32}{=\mathrm{y}}_{32}\right)={\displaystyle \prod _{\mathrm{a}=1}^3{\displaystyle \prod _{\mathrm{b}=1}^3{\mathsf{\pi }}_{\mathrm{ab}}^{{\mathrm{y}}_{\mathrm{ab}}}}}$$

(1)

where ${0<\mathsf{\pi }}_{\mathrm{ab}}<1$; $\mathrm{a},\mathrm{b}=1,2,3$;${\mathrm{y}}_{\mathrm{ab}}=0,1$; ${\displaystyle \sum _{\mathrm{a}=1}^3{\displaystyle \sum _{\mathrm{b}=1}^3{\mathrm{Y}}_{\mathrm{ab}}=1}}$; and ${\displaystyle \sum _{\mathrm{a}=1}^3{\displaystyle \sum _{\mathrm{b}=1}^3{\mathsf{\pi }}_{\mathrm{ab}}=1}}$.${\mathsf{\pi }}_{\mathrm{ab}}=\mathrm{P}\left({\mathrm{Y}}_1{=\mathrm{a},\mathrm{Y}}_2=\mathrm{b}\right)$ is the joint probability of the response variables. ${\mathsf{\pi }}_{\mathrm{a}•}={\displaystyle \sum _{\mathrm{b}=1}^3\mathrm{P}\left({\mathrm{Y}}_1{=\mathrm{a},\mathrm{Y}}_2=\mathrm{b}\right)}$ is the marginal probability for the variable ${\mathrm{Y}}_1$ and ${\mathsf{\pi }}_{•\mathrm{b}}={\displaystyle \sum _{\mathrm{a}=1}^3\mathrm{P}\left({\mathrm{Y}}_1{=\mathrm{a},\mathrm{Y}}_2=\mathrm{b}\right)}$ is the marginal probability for the variable ${\mathrm{Y}}_2$.

The model that can be used for ordinal logistic regression is the cumulative logit model. In this study, the BPOLR model has two ordinal-scale response variables that are correlated with each other and each has 3 categories, so we have the cumulative probability for ${\mathrm{Y}}_1$ and ${\mathrm{Y}}_2,$ i.e., $\mathrm{P}\left({\mathrm{Y}}_1\le \mathrm{a}\left|x\right.\right);\mathrm{a}=1,2$ and $\mathrm{P}\left({\mathrm{Y}}_2\le \mathrm{b}\left|x\right.\right);\mathrm{b}=1,2.$ Let $F_{\mathrm{a}•}\left(x\right)=\mathrm{P}\left({\mathrm{Y}}_1\le \mathrm{a}\left|x\right.\right);\mathrm{a}=1,2$ is the marginal cumulative probability of variable $Y_1$ given $x$, $F_{•\mathrm{b}}\left(x\right)=\mathrm{P}\left({\mathrm{Y}}_2\le \mathrm{b}\left|x\right.\right);\mathrm{b}=1,2$ is the marginal cumulative probability of variable $Y_2$ given $x$ and $F_{\mathrm{ab}}\left(x\right)=\mathrm{P}\left({\mathrm{Y}}_1\le {\mathrm{a},\mathrm{Y}}_2\le \mathrm{b}\left|x\right.\right);\mathrm{a},\mathrm{b}=1,2$ is the joint cumulative probability between variables $Y_1$ and $Y_2$ given $x$, so we have the BPOLR model as follows:

• The cumulative logit model for ${\mathrm{Y}}_1$

  $$g_1=\mathrm{logit}\left(\mathrm{P}\left({\mathrm{Y}}_1\le 1\left|x\right.\right)\right)=\mathrm{logit}\left(F_{1•}\left(x\right)\right)=\ln \left(\frac{F_{1•}\left(x\right)}{1-F_{1•}\left(x\right)}\right)={\alpha }_{11}+{\displaystyle \sum _{j=1}^k{\mathsf{\beta }}_{1j}{}^T}{x}_{j\left(r_j\right)}^{*}$$

  (2)

  $$g_2=\mathrm{logit}\left(\mathrm{P}\left({\mathrm{Y}}_1\le 2\left|x\right.\right)\right)=\mathrm{logit}\left(F_{2•}\left(x\right)\right)=\ln \left(\frac{F_{2•}\left(x\right)}{1-F_{2•}\left(x\right)}\right)={\alpha }_{12}+{\displaystyle \sum _{j=1}^k{\mathsf{\beta }}_{1j}{}^T}{x}_{j\left(r_j\right)}^{*}$$

  (3)
• The cumulative logit model for ${\mathrm{Y}}_2$

  $$g_3=\mathrm{logit}\left(\mathrm{P}\left({\mathrm{Y}}_2\le 1\left|x\right.\right)\right)=\mathrm{logit}\left(F_{•1}\left(x\right)\right)=\ln \left(\frac{F_{•1}\left(x\right)}{1-F_{•1}\left(x\right)}\right)={\alpha }_{21}+{\displaystyle \sum _{j=1}^k{\mathsf{\beta }}_{2j}{}^T}{x}_{j\left(r_j\right)}^{*}$$

  (4)

  $$g_4=\mathrm{logit}\left(\mathrm{P}\left({\mathrm{Y}}_2\le 2\left|x\right.\right)\right)=\mathrm{logit}\left(F_{•2}\left(x\right)\right)=\ln \left(\frac{F_{•2}\left(x\right)}{1-F_{•2}\left(x\right)}\right)={\alpha }_{22}+{\displaystyle \sum _{j=1}^k{\mathsf{\beta }}_{2j}{}^T}{x}_{j\left(r_j\right)}^{*}$$

  (5)
• The odds ratio transformation model for ${\mathrm{Y}}_1$ and ${\mathrm{Y}}_2$

  $$g_5=\ln \left({\psi }_{11}\left(x\right)\right)=\ln \left(\frac{F_{11}\left(x\right)\left[1-F_{1•}\left(x\right)-F_{•1}\left(x\right)+F_{11}\left(x\right)\right]}{\left[F_{1•}\left(x\right)-F_{11}\left(x\right)\right]\left[F_{•1}\left(x\right)-F_{11}\left(x\right)\right]}\right)={\Delta }_{11}+{\displaystyle \sum _{j=1}^k{\mathsf{\gamma }}_j{}^T}{x}_{j\left(r_j\right)}^{*}$$

  (6)

  $$g_6=\ln \left({\psi }_{12}\left(x\right)\right)=\ln \left(\frac{F_{12}\left(x\right)\left[1-F_{1•}\left(x\right)-F_{•2}\left(x\right)+F_{12}\left(x\right)\right]}{\left[F_{1•}\left(x\right)-F_{12}\left(x\right)\right]\left[F_{•2}\left(x\right)-F_{12}\left(x\right)\right]}\right)={\Delta }_{12}+{\displaystyle \sum _{j=1}^k{\mathsf{\gamma }}_j{}^T}{x}_{j\left(r_j\right)}^{*}$$

  (7)

  $$g_7=\ln \left({\psi }_{21}\left(x\right)\right)=\ln \left(\frac{F_{21}\left(x\right)\left[1-F_{2•}\left(x\right)-F_{•1}\left(x\right)+F_{21}\left(x\right)\right]}{\left[F_{2•}\left(x\right)-F_{21}\left(x\right)\right]\left[F_{•1}\left(x\right)-F_{21}\left(x\right)\right]}\right)={\Delta }_{21}+{\displaystyle \sum _{j=1}^k{\mathsf{\gamma }}_j{}^T}{x}_{j\left(r_j\right)}^{*}$$

  (8)

  $$g_8=\ln \left({\psi }_{22}\left(x\right)\right)=\ln \left(\frac{F_{22}\left(x\right)\left[1-F_{2•}\left(x\right)-F_{•2}\left(x\right)+F_{22}\left(x\right)\right]}{\left[F_{2•}\left(x\right)-F_{22}\left(x\right)\right]\left[F_{•2}\left(x\right)-F_{22}\left(x\right)\right]}\right)={\Delta }_{22}+{\displaystyle \sum _{j=1}^k{\mathsf{\gamma }}_j{}^T}{x}_{j\left(r_j\right)}^{*}$$

  (9)

  where $\left\{{\mathsf{\alpha }}_{1\mathrm{a}}{,\mathsf{\alpha }}_{2\mathrm{b}}{,\Delta }_{\mathrm{ab}}\right\} ; \mathrm{a},\mathrm{b}=1,2$ are the intercept parameters with ${\alpha }_{11}<{\alpha }_{12}$ and ${\alpha }_{21}<{\alpha }_{22}$; ${\mathsf{\beta }}_{1j}^{}$,${\mathsf{\beta }}_{2j}^{}$ and ${\mathsf{\gamma }}_j^{}$ are vector of parameters for the j-th predictor variable, which are symbolized by ${\mathsf{\beta }}_{1j}^{}={\left[\begin{array}{cccc}{\beta }_{01j}& {\beta }_{11j}& \dots & {\beta }_{r1j}\end{array}\right]}^T$, ${\mathsf{\beta }}_{2j}^{}={\left[\begin{array}{cccc}{\beta }_{02j}& {\beta }_{12j}& \dots & {\beta }_{r2j}\end{array}\right]}^T$,${\mathsf{\gamma }}_j^{}={\left[\begin{array}{cccc}{\gamma }_{0j}& {\gamma }_{1j}& \dots & {\gamma }_{rj}\end{array}\right]}^T$. $x$ is vector of predictor variable with $x={\left[\begin{array}{cccc}1& x_{1\left(r_1\right)}^{*}& \dots & x_{k\left(r_k\right)}^{*}\end{array}\right]}^T$ where $x_{j\left(r_j\right)}^{*}={\left[\begin{array}{ccccc}1& x_j& x_j^2& \dots & x_j^{r_j}\end{array}\right]}^T$ is vector of the j-th predictor variable with the r-th degree of polynomial.

Based on Equations (2)–(5), the marginal cumulative probability of variable ${\mathrm{Y}}_1$ and ${\mathrm{Y}}_2$ can be obtained from the following equation:

$$\mathrm{P}\left({\mathrm{Y}}_1\le \mathrm{a}\left|x\right.\right)=F_{\mathrm{a}·}\left(x\right)=\frac{\exp \left({\alpha }_{1\mathrm{a}}+{\displaystyle \sum _{j=1}^k{\mathsf{\beta }}_{1j}{}^T}{x}_{j\left(r_j\right)}^{*}\right)}{1+\exp \left({\alpha }_{1\mathrm{a}}+{\displaystyle \sum _{j=1}^k{\mathsf{\beta }}_{1j}{}^T}{x}_{j\left(r_j\right)}^{*}\right)}; \mathrm{a}=1,2$$

(10)

$$\mathrm{P}\left({\mathrm{Y}}_2\le \mathrm{b}\left|x\right.\right)=F_{·\mathrm{b}}\left(x\right)=\frac{\exp \left({\alpha }_{2\mathrm{b}}+{\displaystyle \sum _{j=1}^k{\mathsf{\beta }}_{2j}{}^T}{x}_{j\left(r_j\right)}^{*}\right)}{1+\exp \left({\alpha }_{2\mathrm{b}}+{\displaystyle \sum _{j=1}^k{\mathsf{\beta }}_{2j}{}^T}{x}_{j\left(r_j\right)}^{*}\right)} ; \mathrm{b}=1,2$$

(11)

${\psi }_{\mathrm{ab}}\left(x\right)$ is the odds ratio (OR), a measure of association that can indicate a correlation between the response variables. The OR value is formulated as follows:

$${\psi }_{\mathrm{ab}}\left(x\right)=\frac{\mathrm{P}\left({\mathrm{Y}}_1\le {\mathrm{a},\mathrm{Y}}_2\le \mathrm{b}\left|x\right.\right)\mathrm{P}\left({\mathrm{Y}}_1{>\mathrm{a},\mathrm{Y}}_2>\mathrm{b}\left|x\right.\right)}{\mathrm{P}\left({\mathrm{Y}}_1{>\mathrm{a},\mathrm{Y}}_2\le \mathrm{b}\left|x\right.\right)\mathrm{P}\left({\mathrm{Y}}_1\le {\mathrm{a},\mathrm{Y}}_2>\mathrm{b}\left|x\right.\right)}=\frac{F_{\mathrm{ab}}\left(x\right)\left[1-F_{\mathrm{a}•}\left(x\right)-F_{•\mathrm{b}}\left(x\right)+F_{\mathrm{ab}}\left(x\right)\right]}{\left[F_{\mathrm{a}•}\left(x\right)-F_{\mathrm{ab}}\left(x\right)\right]\left[F_{•\mathrm{b}}\left(x\right)-F_{\mathrm{ab}}\left(x\right)\right]}$$

(12)

where $F_{\mathrm{ab}}\left(x\right)$ can be obtained from:

$$F_{\mathrm{ab}}\left(x\right)=\left\{\begin{array}{cc}\frac{\mathrm{c}-\sqrt{{\mathrm{c}}_{}^2+\mathrm{d}}}{2\left({\psi }_{\mathrm{ab}}\left(x\right)-1\right)}& , {\psi }_{\mathrm{ab}}\left(x\right)\ne 1\\ F_{\mathrm{a}•}\left(x\right)F_{•\mathrm{b}}\left(x\right)& , {\psi }_{\mathrm{ab}}\left(x\right)=1\end{array}\right.$$

(13)

where $\mathrm{a},\mathrm{b}=1,2$; $\mathrm{c}=1+\left({\psi }_{\mathrm{ab}}\left(x\right)-1\right)\left(F_{\mathrm{a}•}\left(\mathrm{x}\right)+F_{•\mathrm{b}}\left(x\right)\right)$ and $\mathrm{d}=-4{\psi }_{\mathrm{ab}}\left(x\right)\left({\psi }_{\mathrm{ab}}\left(x\right)-1\right)F_{\mathrm{a}•}\left(x\right)F_{•\mathrm{b}}\left(x\right)$.

Next, the joint probability of the response variables,${\mathsf{\pi }}_{\mathrm{ab}}$, can be obtained as follows:

${\mathsf{\pi }}_{11}=F_{11}$; ${\mathsf{\pi }}_{12}=F_{12}-F_{11}$; ${\mathsf{\pi }}_{13}=F_{1•}-F_{12}$; ${\mathsf{\pi }}_{21}=F_{21}-F_{11}$; ${\mathsf{\pi }}_{22}=F_{22}-F_{12}-F_{21}+F_{11}$;${\mathsf{\pi }}_{23}=F_{2•}-F_{22}-F_{1•}+F_{12}$; ${\mathsf{\pi }}_{31}=F_{•1}-F_{21}$; ${\mathsf{\pi }}_{32}=F_{•2}-F_{22}-F_{•1}+F_{21}$; ${\mathsf{\pi }}_{33}=1-F_{2•}-F_{•2}+F_{22}$.

3. Parameter Estimation of The BPOLR Model

In this study, the parameter estimation of the BPOLR model was carried out using the maximum likelihood estimation (MLE) method. The principle of the MLE method is to estimate the parameters of the BPOLR model, namely:

$$\begin{gathered} \mathsf{\theta }=\left[\begin{array}{cccccccccccc}{\alpha }_{11}& {\alpha }_{12}& {\beta }_{011}& {\beta }_{111}& \dots & {\beta }_{r1k}& {\alpha }_{21}& {\alpha }_{22}& {\beta }_{021}& {\beta }_{121}& \dots & {\beta }_{r2k}\end{array}\right. \\ {\left.\begin{array}{cccccccc}{\Delta }_{11}& {\Delta }_{12}& {\Delta }_{21}& {\Delta }_{22}& {\gamma }_{01}& {\gamma }_{11}& \dots & {\gamma }_{rk}\end{array}\right]}^T \end{gathered}$$

(14)

obtained by maximizing the likelihood function. Based on Equation (1), the likelihood function is obtained as follows:

$$\mathrm{L}\left(\mathsf{\theta }\right)={\displaystyle \prod _{\mathrm{i}=1}^{\mathrm{n}}{\displaystyle \prod _{\mathrm{a}=1}^3{\displaystyle \prod _{\mathrm{b}=1}^3{\mathsf{\pi }}_{\mathrm{abi}}^{{\mathrm{y}}_{\mathrm{abi}}}}}}$$

(15)

To simplify the calculation, an ln transformation is carried out on the likelihood function so that the ln-likelihood function is formed as follows:

$$\mathrm{Q}=\mathrm{lnL}\left(\mathsf{\theta }\right)={\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{n}}{\displaystyle \sum _{\mathrm{a}=1}^3{\displaystyle \sum _{\mathrm{b}=1}^3{\mathrm{y}}_{\mathrm{abi}}\ln }}}\left({\mathsf{\pi }}_{\mathrm{abi}}^{}\right)={\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{n}}\left[{\mathrm{y}}_{11\mathrm{i}}{\ln \mathsf{\pi }}_{11\mathrm{i}}^{}{+\mathrm{y}}_{12\mathrm{i}}{\ln \mathsf{\pi }}_{12\mathrm{i}}^{}{+\dots +\mathrm{y}}_{33\mathrm{i}}{\ln \mathsf{\pi }}_{33\mathrm{i}}^{}\right]}$$

(16)

the next step of the ln-likelihood function is the first partial derivative of the parameter to be estimated and then equated with zero. The results of the first partial derivative of the ln-likelihood function with respect to its parameters are as follows:

• The first partial derivative of the ln-likelihood function to the parameter ${\alpha }_{11}$

  $$\frac{\partial \mathbb{Q}}{\partial {\alpha }_{11}}={\displaystyle \sum _{i=1}^n\left[\left(\frac{y_{11i}}{{\pi }_{11i}}-\frac{y_{12i}}{{\pi }_{12i}}-\frac{y_{21i}}{{\pi }_{21i}}+\frac{y_{22i}}{{\pi }_{22i}}\right){\delta }_{11i}+\left(\frac{y_{12i}}{{\pi }_{12i}}-\frac{y_{13i}}{{\pi }_{13i}}-\frac{y_{22i}}{{\pi }_{22i}}+\frac{y_{23i}}{{\pi }_{23i}}\right){\delta }_{12i}+\left(\frac{y_{13i}}{{\pi }_{13i}}-\frac{y_{23i}}{{\pi }_{23i}}\right)\right]}{\xi }_{1•i}$$

  (17)
• The first partial derivative of the ln-likelihood function to the parameter ${\alpha }_{12}$

  $$\frac{\partial \mathbb{Q}}{\partial {\alpha }_{12}}={\displaystyle \sum _{i=1}^n\left[\left(\frac{y_{21i}}{{\pi }_{21i}}-\frac{y_{22i}}{{\pi }_{22i}}-\frac{y_{31i}}{{\pi }_{31i}}+\frac{y_{32i}}{{\pi }_{32i}}\right){\delta }_{21i}+\left(\frac{y_{22i}}{{\pi }_{22i}}-\frac{y_{23i}}{{\pi }_{23i}}-\frac{y_{32i}}{{\pi }_{32i}}+\frac{y_{33i}}{{\pi }_{33i}}\right){\delta }_{22i}+\left(\frac{y_{23i}}{{\pi }_{23i}}-\frac{y_{33i}}{{\pi }_{33i}}\right)\right]}{\xi }_{2•i}$$

  (18)
• The first partial derivative of the ln-likelihood function to the parameter ${\mathsf{\beta }}_{1j}; j=1,2,\dots ,k$

  $$\begin{gathered} \frac{\partial \mathbb{Q}}{\partial {\mathsf{\beta }}_{1j}{}^{}}={\displaystyle \sum _{i=1}^n\left[\left\{\left(\frac{y_{11i}}{{\pi }_{11i}}-\frac{y_{12i}}{{\pi }_{12i}}-\frac{y_{21i}}{{\pi }_{21i}}+\frac{y_{22i}}{{\pi }_{22i}}\right){\delta }_{11i}+\left(\frac{y_{12i}}{{\pi }_{12i}}-\frac{y_{13i}}{{\pi }_{13i}}-\frac{y_{22i}}{{\pi }_{22i}}+\frac{y_{23i}}{{\pi }_{23i}}\right){\delta }_{12i}+\left(\frac{y_{13i}}{{\pi }_{13i}}-\frac{y_{23i}}{{\pi }_{23i}}\right)\right\}{\xi }_{1•i}\right.}{x}_{j\left(r_j\right)}^{*}+ \\ \left.\left\{\left(\frac{y_{21i}}{{\pi }_{21i}}-\frac{y_{22i}}{{\pi }_{22i}}-\frac{y_{31i}}{{\pi }_{31i}}+\frac{y_{32i}}{{\pi }_{32i}}\right){\delta }_{21i}+\left(\frac{y_{22i}}{{\pi }_{22i}}-\frac{y_{23i}}{{\pi }_{23i}}-\frac{y_{32i}}{{\pi }_{32i}}+\frac{y_{33i}}{{\pi }_{33i}}\right){\delta }_{22i}+\left(\frac{y_{23i}}{{\pi }_{23i}}-\frac{y_{33i}}{{\pi }_{33i}}\right)\right\}{\xi }_{2•i}{x}_{j\left(r_j\right)}^{*}\right] \end{gathered}$$

  (19)
• The first partial derivative of the ln-likelihood function to the parameter ${\alpha }_{21}$

  $$\frac{\partial \mathbb{Q}}{\partial {\alpha }_{21}}={\displaystyle \sum _{i=1}^n\left[\left(\frac{y_{11i}}{{\pi }_{11i}}-\frac{y_{12i}}{{\pi }_{12i}}-\frac{y_{21i}}{{\pi }_{21i}}+\frac{y_{22i}}{{\pi }_{22i}}\right){\mu }_{11i}+\left(\frac{y_{21i}}{{\pi }_{21i}}-\frac{y_{22i}}{{\pi }_{22i}}-\frac{y_{31i}}{{\pi }_{31i}}+\frac{y_{32i}}{{\pi }_{32i}}\right){\mu }_{21i}+\left(\frac{y_{31i}}{{\pi }_{31i}}-\frac{y_{32i}}{{\pi }_{32i}}\right)\right]}{\xi }_{•1i}$$

  (20)
• The first partial derivative of the ln-likelihood function to the parameter ${\alpha }_{22}$

  $$\frac{\partial \mathbb{Q}}{\partial {\alpha }_{22}}={\displaystyle \sum _{i=1}^n\left[\left(\frac{y_{12i}}{{\pi }_{12i}}-\frac{y_{13i}}{{\pi }_{13i}}-\frac{y_{22i}}{{\pi }_{22i}}+\frac{y_{23i}}{{\pi }_{23i}}\right){\mu }_{12i}+\left(\frac{y_{22i}}{{\pi }_{22i}}-\frac{y_{23i}}{{\pi }_{23i}}-\frac{y_{32i}}{{\pi }_{32i}}+\frac{y_{33i}}{{\pi }_{33i}}\right){\mu }_{22i}+\left(\frac{y_{32i}}{{\pi }_{32i}}-\frac{y_{33i}}{{\pi }_{33i}}\right)\right]}{\xi }_{•2i}$$

  (21)
• The first partial derivative of the ln-likelihood function to the parameter ${\mathsf{\beta }}_{2j}; j=1,2,\dots ,k$

  $$\begin{gathered} \frac{\partial \mathbb{Q}}{\partial {\mathsf{\beta }}_{2j}}={\displaystyle \sum _{i=1}^n\left[\left\{\left(\frac{y_{11i}}{{\pi }_{11i}}-\frac{y_{12i}}{{\pi }_{12i}}-\frac{y_{21i}}{{\pi }_{21i}}+\frac{y_{22i}}{{\pi }_{22i}}\right){\mu }_{11i}+\left(\frac{y_{21i}}{{\pi }_{21i}}-\frac{y_{22i}}{{\pi }_{22i}}-\frac{y_{31i}}{{\pi }_{31i}}+\frac{y_{32i}}{{\pi }_{32i}}\right){\mu }_{21i}+\left(\frac{y_{31i}}{{\pi }_{31i}}-\frac{y_{32i}}{{\pi }_{32i}}\right)\right\}{\xi }_{•1i}\right.}{x}_{j\left(r_j\right)}^{*}+ \\ \left.\left\{\left(\frac{y_{12i}}{{\pi }_{12i}}-\frac{y_{13i}}{{\pi }_{13i}}-\frac{y_{22i}}{{\pi }_{22i}}+\frac{y_{23i}}{{\pi }_{23i}}\right){\mu }_{12i}+\left(\frac{y_{22i}}{{\pi }_{22i}}-\frac{y_{23i}}{{\pi }_{23i}}-\frac{y_{32i}}{{\pi }_{32i}}+\frac{y_{33i}}{{\pi }_{33i}}\right){\mu }_{22i}+\left(\frac{y_{32i}}{{\pi }_{32i}}-\frac{y_{33i}}{{\pi }_{33i}}\right)\right\}{\xi }_{•2i}{x}_{j\left(r_j\right)}^{*}\right] \end{gathered}$$

  (22)
• The first partial derivative of the ln-likelihood function to the parameter ${\Delta }_{11}$

  $$\frac{\partial \mathbb{Q}}{\partial {\Delta }_{11}}={\displaystyle \sum _{i=1}^n\left[\left(\frac{y_{11i}}{{\pi }_{11i}}-\frac{y_{12i}}{{\pi }_{12i}}-\frac{y_{21i}}{{\pi }_{21i}}+\frac{y_{22i}}{{\pi }_{22i}}\right)z_{11i}{\psi }_{11i}\right]}$$

  (23)
• The first partial derivative of the ln-likelihood function to the parameter ${\Delta }_{12}$

  $$\frac{\partial \mathbb{Q}}{\partial {\Delta }_{12}}={\displaystyle \sum _{i=1}^n\left[\left(\frac{y_{12i}}{{\pi }_{12i}}-\frac{y_{13i}}{{\pi }_{13i}}-\frac{y_{22i}}{{\pi }_{22i}}+\frac{y_{23i}}{{\pi }_{23i}}\right)z_{12i}{\psi }_{12i}\right]}$$

  (24)
• The first partial derivative of the ln-likelihood function to the parameter ${\Delta }_{21}$

  $$\frac{\partial \mathbb{Q}}{\partial {\Delta }_{21}}={\displaystyle \sum _{i=1}^n\left[\left(\frac{y_{21i}}{{\pi }_{21i}}-\frac{y_{22i}}{{\pi }_{22i}}-\frac{y_{31i}}{{\pi }_{31i}}+\frac{y_{32i}}{{\pi }_{32i}}\right)z_{21i}{\psi }_{21i}\right]}$$

  (25)
• The first partial derivative of the ln-likelihood function to the parameter ${\Delta }_{22}$

  $$\frac{\partial \mathbb{Q}}{\partial {\Delta }_{22}}={\displaystyle \sum _{i=1}^n\left[\left(\frac{y_{22i}}{{\pi }_{22i}}-\frac{y_{23i}}{{\pi }_{23i}}-\frac{y_{32i}}{{\pi }_{32i}}+\frac{y_{33i}}{{\pi }_{33i}}\right)z_{22i}{\psi }_{22i}\right]}$$

  (26)
• The first partial derivative of the ln-likelihood function to the parameter ${\mathsf{\gamma }}_j; j=1,2,\dots ,k$

  $$\begin{gathered} \frac{\partial \mathbb{Q}}{\partial {\mathsf{\gamma }}_j}={\displaystyle \sum _{i=1}^n\left[\left\{\left(\frac{y_{11i}}{{\pi }_{11i}}-\frac{y_{12i}}{{\pi }_{12i}}-\frac{y_{21i}}{{\pi }_{21i}}+\frac{y_{22i}}{{\pi }_{22i}}\right)z_{11i}{\psi }_{11i}+\left(\frac{y_{12i}}{{\pi }_{12i}}-\frac{y_{13i}}{{\pi }_{13i}}-\frac{y_{22i}}{{\pi }_{22i}}+\frac{y_{23i}}{{\pi }_{23i}}\right)z_{12i}{\psi }_{12i}+\right.\right.} \\ \left.\left.\left(\frac{y_{21i}}{{\pi }_{21i}}-\frac{y_{22i}}{{\pi }_{22i}}-\frac{y_{31i}}{{\pi }_{31i}}+\frac{y_{32i}}{{\pi }_{32i}}\right)z_{21i}{\psi }_{21i}+\left(\frac{y_{22i}}{{\pi }_{22i}}-\frac{y_{23i}}{{\pi }_{23i}}-\frac{y_{32i}}{{\pi }_{32i}}+\frac{y_{33i}}{{\pi }_{33i}}\right)z_{22i}{\psi }_{22i}\right\}{x}_{j\left(r_j\right)}^{*}\right] \end{gathered}$$

  (27)

  where ${\delta }_{\mathrm{ab}}=\frac{1}{2}\left\{1-{\mathrm{S}}_{\mathrm{ab}}{}^{-1}\left[1+{\psi }_{\mathrm{ab}}\left(F_{\mathrm{a}•}-F_{•\mathrm{b}}\right)-F_{\mathrm{a}•}-F_{•\mathrm{b}}\right]\right\}$; ${\mu }_{\mathrm{ab}}=\frac{1}{2}\left\{1-{\mathrm{S}}_{\mathrm{ab}}{}^{-1}\left[1+{\psi }_{\mathrm{ab}}\left(F_{•\mathrm{b}}-F_{\mathrm{a}•}\right)-F_{\mathrm{a}•}-F_{•\mathrm{b}}\right]\right\}$

${\mathrm{S}}_{\mathrm{ab}}={\left[{\left[1+\left({\psi }_{\mathrm{ab}}-1\right)\left(F_{\mathrm{a}•}+F_{•\mathrm{b}}\right)\right]}^2-4{\psi }_{\mathrm{ab}}\left({\psi }_{\mathrm{ab}}-1\right)F_{\mathrm{a}•}{F}_{•\mathrm{b}}\right]}^{\frac{1}{2}}$; ${\xi }_{\mathrm{a}•}=F_{\mathrm{a}•}\left(1-F_{\mathrm{a}•}\right)$; ${\xi }_{•\mathrm{b}}=F_{•\mathrm{b}}\left(1-F_{•\mathrm{b}}\right)$

$z_{\mathrm{ab}}={\left[2{\left({\psi }_{\mathrm{ab}}-1\right)}^2{\mathrm{S}}_{\mathrm{ab}}\right]}^{-1}\left\{1-{\mathrm{S}}_{\mathrm{ab}}+\left({\psi }_{\mathrm{ab}}-1\right)\left[F_{\mathrm{a}•}+F_{•\mathrm{b}}-2F_{\mathrm{a}•}{F}_{•\mathrm{b}}\right]\right\}$.

The result of the first partial derivatives are not closed-form, thus, a numerical optimization such as the Berndt–Hall–Hall–Hausman (BHHH) method is needed to obtain the maximum likelihood estimators and the BHHH algorithm is as follows:

• Step 1. Determine the initial value for

  $${\mathsf{\theta }}^{(0)}=\left[\begin{array}{cccccccccccc}{\alpha }_{11}^{(0)}& {\alpha }_{12}^{(0)}& {\beta }_{011}^{(0)}& {\beta }_{111}^{(0)}& \dots & {\beta }_{r1k}^{(0)}& {\alpha }_{21}^{(0)}& {\alpha }_{22}^{(0)}& {\beta }_{021}^{(0)}& {\beta }_{121}^{(0)}& \dots & {\beta }_{r2k}^{(0)}\end{array}\right.$$

  $${\left.\begin{array}{cccccccc}{\Delta }_{11}^{(0)}& {\Delta }_{12}^{(0)}& {\Delta }_{21}^{(0)}& {\Delta }_{22}^{(0)}& {\gamma }_{01}^{(0)}& {\gamma }_{11}^{(0)}& \dots & {\gamma }_{rk}^{(0)}\end{array}\right]}^T$$

  obtained from the parameter estimator of the ordinal logistic regression model on each response variable.
• Step 2. Calculate the gradient vector elements obtained from the first partial derivative of the ln-likelihood function for each parameter

  $$q\left(\mathsf{\theta }\right)=\left[\begin{array}{cccccccccccc}\frac{\partial \mathbb{Q}}{\partial {\alpha }_{11}}& \frac{\partial \mathbb{Q}}{\partial {\alpha }_{12}}& \frac{\partial \mathbb{Q}}{\partial {\beta }_{011}}& \frac{\partial \mathbb{Q}}{\partial {\beta }_{111}}& \dots & \frac{\partial \mathbb{Q}}{\partial {\beta }_{r1k}}& \frac{\partial \mathbb{Q}}{\partial {\alpha }_{21}}& \frac{\partial \mathbb{Q}}{\partial {\alpha }_{22}}& \frac{\partial \mathbb{Q}}{\partial {\beta }_{021}}& \frac{\partial \mathbb{Q}}{\partial {\beta }_{121}}& \dots & \frac{\partial \mathbb{Q}}{\partial {\beta }_{r2k}}\end{array}\right.$$

  $${\left.\begin{array}{cccccccc}\frac{\partial \mathbb{Q}}{\partial {\Delta }_{11}}& \frac{\partial \mathbb{Q}}{\partial {\Delta }_{12}}& \frac{\partial \mathbb{Q}}{\partial {\Delta }_{21}}& \frac{\partial \mathbb{Q}}{\partial {\Delta }_{22}}& \frac{\partial \mathbb{Q}}{\partial {\gamma }_{01}}& \frac{\partial \mathbb{Q}}{\partial {\gamma }_{11}}& \dots & \frac{\partial \mathbb{Q}}{\partial {\gamma }_{rk}}\end{array}\right]}^T$$
• Step 3. Calculate the Hessian matrix $H\left(\mathsf{\theta }\right)$ that can be obtained from the following formula

  $$H\left(\mathsf{\theta }\right)=-\left[q^T\left(\mathsf{\theta }\right)q\left(\mathsf{\theta }\right)\right]$$
• Step 4. Start the BHHH iteration process with the following formula

  $${\widehat{\mathsf{\theta }}}^{(t+1)}={\widehat{\mathsf{\theta }}}^{(t)}-H^{-1}\left({\widehat{\mathsf{\theta }}}^{(t)}\right)q\left({\widehat{\mathsf{\theta }}}^{(t)}\right), \mathrm{for} t=0,1,2,....$$
• Step 5. The iteration will stop if $‖{\widehat{\mathsf{\theta }}}_{}^{(t+1)}-{\widehat{\mathsf{\theta }}}_{}^{(t)}‖<\mathsf{\epsilon }$, where $\mathsf{\epsilon }$ is a very small positive number. The last iteration produces an estimator value for each parameter.

Furthermore, the best model in this study is determined using the criteria of Akaike information criterion corrected (AICc), Bayesian information criterion (BIC), and deviance which are defined as follows:

$$AICc=-2\ln L\left(\widehat{\mathsf{\theta }}\right)+2K+\frac{2K\left(K+1\right)}{n-K-1}$$

(28)

$$BIC=-2\ln L\left(\widehat{\mathsf{\theta }}\right)+K\ln n$$

(29)

$$Deviance=-2\ln L\left(\widehat{\mathsf{\theta }}\right)$$

(30)

where $L\left(\widehat{\mathsf{\theta }}\right)$ is the likelihood value of the parameter’s estimate,$K$ is the number of covariates, and 𝑛 is the sample size. The best model is the BPOLR model, which has the smallest values of AICc, BIC, and deviance.

4. Hypothesis Testing of The BPOLR Model

The simultaneous hypothesis testing of the BPOLR model aims to determine the influence of the predictor variable on the response variable simultaneously or to discover at least one predictor variable which has a significant effect against the response variable. The hypothesis used is as follows:

$$H_0:{\beta }_{p1q}={\beta }_{p2q}={\gamma }_{pq}=0 ; p=1,2,\dots ,r_q ; q=1,2,\dots ,k$$

(31)

$$H_1: \mathrm{at} \mathrm{least} \mathrm{one} {\beta }_{p1q}\ne 0 \mathrm{or} {\beta }_{p2q}\ne 0 \mathrm{or} {\gamma }_{pq}\ne 0$$

In this study, we used the MLRT method to obtain the test statistics by calculating the ratio between the maximum value of the likelihood function under $H_0$, $L\left(\omega \right)$, with the maximum value of the likelihood function under population, $L\left(\mathsf{\Omega }\right)$.

The first step to obtaining the maximum value of $L\left(\omega \right)$ is to define the parameter space under $H_0$, $\omega =\left\{{\alpha }_{11},{\alpha }_{12},{\alpha }_{21},{\alpha }_{22},{\Delta }_{11},{\Delta }_{12},{\Delta }_{21},{\Delta }_{22},{\beta }_{011},\dots ,{\beta }_{01k},{\beta }_{021},\dots ,{\beta }_{02k},\right.$$\left.{\gamma }_{01},\dots ,{\gamma }_{0k}\right\}$. The likelihood function under $H_0$ is:

$$L\left(\omega \right)={\displaystyle \prod _{i=1}^n{\displaystyle \prod _{a=1}^3{\displaystyle \prod _{b=1}^3{\left({\pi }_{abi}^{*}\right)}^{y_{abi}}}}}$$

(32)

Let ${\mathbb{Q}}_{\omega }=\ln L\left(\omega \right)$, then the log-likelihood function under $H_0$ is

$${\mathbb{Q}}_{\omega }={\displaystyle \sum _{i=1}^n\left(y_{11i}\ln {\pi }_{11i}^{*}+y_{12i}\ln {\pi }_{12i}^{*}+\dots +y_{33i}\ln {\pi }_{33i}^{*}\right)}$$

(33)

Furthermore, we can obtain the estimated value of the maximum log-likelihood function under $H_0$:

$${\mathbb{Q}}_{\widehat{\omega }}={\displaystyle \sum _{i=1}^n\left(y_{11i}\ln {\widehat{\pi }}_{11i}^{*}+y_{12i}\ln {\widehat{\pi }}_{12i}^{*}+\dots +y_{33i}\ln {\widehat{\pi }}_{33i}^{*}\right)}$$

(34)

where

$${\widehat{\pi }}_{11i}^{*}={\widehat{F}}_{11i}^{*} ; {\widehat{\pi }}_{12i}^{*}={\widehat{F}}_{12i}^{*}-{\widehat{F}}_{11i}^{*}; {\widehat{\pi }}_{13i}^{*}={\widehat{F}}_{1•i}^{*}-{\widehat{F}}_{12i}^{*}; {\widehat{\pi }}_{21i}^{*}={\widehat{F}}_{21i}^{*}-{\widehat{F}}_{11i}^{*};$$

$${\widehat{\pi }}_{22i}^{*}={\widehat{F}}_{22i}^{*}-{\widehat{F}}_{12i}^{*}-{\widehat{F}}_{21i}^{*}+{\widehat{F}}_{11i}^{*} ;{\widehat{\pi }}_{23i}^{*}={\widehat{F}}_{2•i}^{*}-{\widehat{F}}_{1•i}^{*}-{\widehat{F}}_{22i}^{*}+{\widehat{F}}_{12i}^{*};{\widehat{\pi }}_{31i}^{*}={\widehat{F}}_{•1i}^{*}-{\widehat{F}}_{21i}^{*};$$

$${\widehat{\pi }}_{32i}^{*}={\widehat{F}}_{•2i}^{*}-{\widehat{F}}_{22i}^{*}-{\widehat{F}}_{•1i}^{*}+{\widehat{F}}_{21i}^{*}; {\widehat{\pi }}_{33i}^{*}=1-{\widehat{F}}_{2•i}^{*}-{\widehat{F}}_{•2i}^{*}+{\widehat{F}}_{22i}^{*}$$

with

$${\hat{\mathrm{F}}}_{\mathrm{abi}}^{*}=\frac{{\mathrm{p}}_3-\sqrt{{\mathrm{p}}_3{}_{}^2{+\mathrm{q}}_3}}{2\left({\hat{\mathsf{\psi }}}_{\mathrm{ab}}-1\right)},{\hat{\mathsf{\psi }}}_{\mathrm{ab}}\ne 1; {\mathrm{p}}_3=1+\left({\hat{\mathsf{\psi }}}_{\mathrm{ab}}-1\right)\left({\hat{\mathrm{F}}}_{\mathrm{a}•\mathrm{i}}^{*}{+\hat{\mathrm{F}}}_{•\mathrm{bi}}^{*}\right); {\mathrm{q}}_3=-4{\hat{\mathsf{\psi }}}_{\mathrm{ab}}\left({\hat{\mathsf{\psi }}}_{\mathrm{ab}}-1\right){\hat{\mathrm{F}}}_{\mathrm{a}•\mathrm{i}}^{*}{\hat{\mathrm{F}}}_{•\mathrm{bi}}^{*}$$

$${\hat{\mathrm{F}}}_{\mathrm{a}•\mathrm{i}}^{*}=\frac{\exp \left({\hat{\mathsf{\alpha }}}_{1\mathrm{a}}+{\displaystyle \sum _{\mathrm{j}=1}^{\mathrm{k}}{\hat{\mathsf{\beta }}}_{01\mathrm{j}}}\right)}{\exp \left({\hat{\mathsf{\alpha }}}_{1\mathrm{a}}+{\displaystyle \sum _{\mathrm{j}=1}^{\mathrm{k}}{\hat{\mathsf{\beta }}}_{01\mathrm{j}}}\right)};{\hat{\mathrm{F}}}_{•\mathrm{bi}}^{*}=\frac{\exp \left({\hat{\mathsf{\alpha }}}_{2\mathrm{b}}+{\displaystyle \sum _{\mathrm{j}=1}^{\mathrm{k}}{\hat{\mathsf{\beta }}}_{02\mathrm{j}}}\right)}{\exp \left({\hat{\mathsf{\alpha }}}_{2\mathrm{b}}+{\displaystyle \sum _{\mathrm{j}=1}^{\mathrm{k}}{\hat{\mathsf{\beta }}}_{02\mathrm{j}}}\right)};{\hat{\mathsf{\psi }}}_{11}=\frac{{\mathrm{n}}_{11\mathrm{i}}\left({\mathrm{n}}_{22\mathrm{i}}{+\mathrm{n}}_{23\mathrm{i}}{+\mathrm{n}}_{32\mathrm{i}}{+\mathrm{n}}_{33\mathrm{i}}\right)}{\left({\mathrm{n}}_{12\mathrm{i}}{+\mathrm{n}}_{13\mathrm{i}}\right)\left({\mathrm{n}}_{21\mathrm{i}}{+\mathrm{n}}_{31\mathrm{i}}\right)} ;$$

$${\hat{\mathsf{\psi }}}_{12}=\frac{\left({\mathrm{n}}_{11\mathrm{i}}{+\mathrm{n}}_{12\mathrm{i}}\right)\left({\mathrm{n}}_{23\mathrm{i}}{+\mathrm{n}}_{33\mathrm{i}}\right)}{{\mathrm{n}}_{13\mathrm{i}}\left({\mathrm{n}}_{21\mathrm{i}}{+\mathrm{n}}_{22\mathrm{i}}{+\mathrm{n}}_{31\mathrm{i}}{+\mathrm{n}}_{32\mathrm{i}}\right)};{\hat{\mathsf{\psi }}}_{21}=\frac{\left({\mathrm{n}}_{11\mathrm{i}}{+\mathrm{n}}_{21\mathrm{i}}\right)\left({\mathrm{n}}_{32\mathrm{i}}{+\mathrm{n}}_{33\mathrm{i}}\right)}{\left({\mathrm{n}}_{12\mathrm{i}}{+\mathrm{n}}_{13\mathrm{i}}{+\mathrm{n}}_{22\mathrm{i}}{+\mathrm{n}}_{23\mathrm{i}}\right){\mathrm{n}}_{31\mathrm{i}}};{\hat{\mathsf{\psi }}}_{22}=\frac{\left({\mathrm{n}}_{11\mathrm{i}}{+\mathrm{n}}_{12\mathrm{i}}{+\mathrm{n}}_{21\mathrm{i}}{+\mathrm{n}}_{22\mathrm{i}}\right){\mathrm{n}}_{33\mathrm{i}}}{\left({\mathrm{n}}_{13\mathrm{i}}{+\mathrm{n}}_{23\mathrm{i}}\right)\left({\mathrm{n}}_{31\mathrm{i}}{+\mathrm{n}}_{32\mathrm{i}}\right)}.$$

Furthermore, to obtain the maximum value of $L\left(\mathsf{\Omega }\right)$ is to define the parameter space under population, $\mathsf{\Omega }=\left\{{\alpha }_{11},{\alpha }_{12},{\alpha }_{21},{\alpha }_{22},{\Delta }_{11},{\Delta }_{12},{\Delta }_{21},{\Delta }_{22},{\beta }_{p1q},{\beta }_{p2q},{\gamma }_{pq};\right.p=0,1,2,\dots ,r_q;$$\left.q=1,2,\dots ,k\right\}$. Then, the likelihood function under population is

$$L\left(\mathsf{\Omega }\right)={\displaystyle \prod _{i=1}^n{\displaystyle \prod _{a=1}^3{\displaystyle \prod _{b=1}^3{\left({\pi }_{abi}^{}\right)}^{y_{abi}}}}}$$

(35)

Let ${\mathbb{Q}}_{\mathsf{\Omega }}=\ln L\left(\mathsf{\Omega }\right)$, then the log-likelihood function under population is

$${\mathbb{Q}}_{\mathsf{\Omega }}={\displaystyle \sum _{i=1}^n\left(y_{11i}\ln {\pi }_{11i}^{}+y_{12i}\ln {\pi }_{12i}^{}+\dots +y_{33i}\ln {\pi }_{33i}^{}\right)}$$

(36)

Furthermore, we can obtain the estimated value of the maximum log-likelihood function under population, that is:

$${\mathbb{Q}}_{\widehat{\mathsf{\Omega }}}={\displaystyle \sum _{i=1}^n\left(y_{11i}\ln {\widehat{\pi }}_{11i}^{}+y_{12i}\ln {\widehat{\pi }}_{12i}^{}+\dots +y_{33i}\ln {\widehat{\pi }}_{33i}^{}\right)}$$

(37)

The test statistic using MLRT method is obtained by

$$G^2=-2\ln \left(\frac{L\left(\widehat{\omega }\right)}{L\left(\widehat{\mathsf{\Omega }}\right)}\right)=2\left({\mathbb{Q}}_{\widehat{\mathsf{\Omega }}}-{\mathbb{Q}}_{\widehat{\omega }}\right)$$

(38)

where ${\mathbb{Q}}_{\widehat{\omega }}$ and ${\mathbb{Q}}_{\widehat{\mathsf{\Omega }}}$ are given in Equations (36) and (39), respectively. Furthermore, the test statistic in Equation (38) can be expressed as follows

$$\begin{array}{l}{G}^2=2\left[{\displaystyle \sum _{i=1}^n\left(y_{11i}\ln \left(\frac{{\widehat{\pi }}_{11i}^{}}{{\widehat{\pi }}_{11i}^{*}}\right)+y_{12i}\ln \left(\frac{{\widehat{\pi }}_{12i}^{}}{{\widehat{\pi }}_{12i}^{*}}\right)+\dots +y_{33i}\ln \left(\frac{{\widehat{\pi }}_{33i}^{}}{{\widehat{\pi }}_{33i}^{*}}\right)\right)}\right]\\ =2\left[{\displaystyle \sum _{i=1}^n\left(n{\widehat{\pi }}_{11i}^{}\ln \left(1+\frac{{\widehat{\pi }}_{11i}^{}-{\widehat{\pi }}_{11i}^{*}}{{\widehat{\pi }}_{11i}^{*}}\right)\right)}+{\displaystyle \sum _{i=1}^n\left(n{\widehat{\pi }}_{12i}^{}\ln \left(1+\frac{{\widehat{\pi }}_{12i}^{}-{\widehat{\pi }}_{12i}^{*}}{{\widehat{\pi }}_{12i}^{*}}\right)\right)}+\dots +{\displaystyle \sum _{i=1}^n\left(n{\widehat{\pi }}_{33i}^{}\ln \left(1+\frac{{\widehat{\pi }}_{33i}^{}-{\widehat{\pi }}_{33i}^{*}}{{\widehat{\pi }}_{33i}^{*}}\right)\right)}\right]\end{array}$$

Based on the Taylor expansion of $\ln \left(1+\frac{{\widehat{\pi }}_{abi}^{}-{\widehat{\pi }}_{abi}^{*}}{{\widehat{\pi }}_{abi}^{*}}\right)$ with $a,b=1,2,3$ and $i=1,2,\dots ,n$ it will be proved that $\underset{n\to \infty }{\lim }{G}^2=G_{}^{*}$.

$$\begin{array}{l}{G}_{}^{*}={\displaystyle \sum _{i=1}^n\frac{{\left(n{\widehat{\pi }}_{11i}^{}-n{\widehat{\pi }}_{11i}^{*}\right)}^2}{n{\widehat{\pi }}_{11i}^{*}}}+{\displaystyle \sum _{i=1}^n\frac{{\left(n{\widehat{\pi }}_{12i}^{}-n{\widehat{\pi }}_{12i}^{*}\right)}^2}{n{\widehat{\pi }}_{12i}^{*}}}+\dots +{\displaystyle \sum _{i=1}^n\frac{{\left(n{\widehat{\pi }}_{33i}^{}-n{\widehat{\pi }}_{33i}^{*}\right)}^2}{n{\widehat{\pi }}_{33i}^{*}}}\\ ={\displaystyle \sum _{i=1}^n{\displaystyle \sum _{a=1}^3{\displaystyle \sum _{b=1}^3\frac{{\left(n{\widehat{\pi }}_{abi}^{}-n{\widehat{\pi }}_{abi}^{*}\right)}^2}{n{\widehat{\pi }}_{abi}^{*}}}}}\\ ={\displaystyle \sum _{i=1}^n{\displaystyle \sum _{a=1}^3{\displaystyle \sum _{b=1}^3\frac{{\left(y_{abi}^{}-n{\widehat{\pi }}_{abi}^{*}\right)}^2}{n{\widehat{\pi }}_{abi}^{*}}}}}\end{array}$$

(39)

The test statistic in Equation (39) has an asymptotic Chi-square limit distribution with v degree of freedom, where v is the difference between the number of parameters in the BPOLR model under population and $H_0$, $v=\left(8+3\left({\displaystyle \sum _{j=1}^k{r}_j}+k\right)\right)-\left(8+3k\right)=3{\displaystyle \sum _{j=1}^k{r}_j}$. Therefore, reject $H_0$ when $G^2>{\chi }_{\left(v,1-\alpha \right)}^2$, where ${\chi }_{\left(v,1-\alpha \right)}^2$ is the $1-\alpha$ quantile from a Chi-square distribution (${\chi }^2$) with 𝑣 degree of freedom.

5. Simulation Study

In this section, we outline a simulation study for the present performance of parameter estimation of the BPOLR model based on the BHHH algorithm. The simulation was constructed by two ordinal response variables that each have three categories and three predictor variables. In this paper, the simulation is carried out for one example of the BPOLR model, i.e., the BPOLR model with degree polynomial (1,2,1). The simulation was carried out for 100 replications with different sample sizes, i.e., n = 100, n = 200, and n = 300. The simulation scenarios are decided as follows:

• Generate three predictor variables (X₁, X₂ and X₃) that are constructed from a standard uniform distribution
• Set the initial coefficients of the BPOLR model as follows:

  $${\alpha }_{11}=0.05 ; {\alpha }_{12}=0.5 ; {\beta }_{111}=-1 ; {\beta }_{112}=0.5 ; {\beta }_{212}=2 ; {\beta }_{113}=1.5 ; {\alpha }_{21}=0.25 ; {\alpha }_{22}=1.25 ;$$

  $${\beta }_{121}=-2 ; {\beta }_{122}=2 ; {\beta }_{222}=2 ; {\beta }_{123}=2 ; {\Delta }_{11}=0.1 ; {\Delta }_{12}=0.15 ; {\Delta }_{21}=0.25 ; {\Delta }_{22}=0.3 ;$$

  $${\gamma }_{11}=-1.5 ; {\gamma }_{12}=2 ; {\gamma }_{22}=3 ; {\gamma }_{13}=2$$
• Generate two ordinal response variables (Y₁ and Y₂) with the following steps:

  ▪

  Determine the cumulative logit model for Y₁ and Y₂ as in Equations (2)–(5) and the odds ratio transformation model, as in Equations (6)–(9)

  ▪

  Determine the marginal cumulative probability for Y₁ and Y₂ as in Equations (10) and (11) and the joint cumulative probability as in Equation (13).

  ▪

  Determine the joint probability of Y₁ and Y₂

  ▪

  Generate two ordinal response variables based on the joint probabilty obtained

• Examine the independence of the response variables using the Mantel–Haenszel test to fulfill the assumption of dependence between the response variables in the bivariate model. If the response variable is independent, then the data generation process is repeated until the dependent response variable is obtained.
• Estimate the parameters of the BPOLR model based on the BHHH algorithm
• Repeat the process for up to 100 replications for each sample size
• Calculate the mean of parameter estimated and its standard error (SE)

The results and comments regarding the simulation study are given in the following tables and figures below. The mean of the estimated parameters and their corresponding standard error (SE) for different sample sizes are presented in

Table 2. Mean of the estimated parameter and the corresponding standard errors of the BPOLR model for different sample sizes.

No.	Parameter	True Coeff.	Mean of the est. Parameter			Standard Error
			n = 100	n = 200	n = 300	n = 100	n = 200	n = 300
1	${\alpha }_{11}$	0.05	0.138	0.097	0.005	1.193	0.662	0.519
2	${\alpha }_{12}$	0.5	0.615	0.549	0.474	1.195	0.666	0.521
3	${\beta }_{111}$	−1	−1.069	−1.176	−0.96	1.144	0.632	0.501
4	${\beta }_{112}$	0.5	−0.005	0.679	0.283	4.856	2.566	2.004
5	${\beta }_{212}$	2	3.213	1.949	2.418	5.492	2.762	2.169
6	${\beta }_{113}$	1.5	1.607	1.604	1.582	1.159	0.652	0.516
7	${\alpha }_{21}$	0.25	0.382	0.337	0.236	1.556	0.777	0.604
8	${\alpha }_{22}$	1.25	1.35	1.349	1.282	1.573	0.789	0.614
9	${\beta }_{121}$	−2	−2.34	−2.216	−2.032	1.546	0.826	0.625
10	${\beta }_{122}$	2	1.119	1.686	1.873	6.949	3.417	2.513
11	${\beta }_{222}$	2	4.466	2.928	2.47	9.294	4.291	3
12	${\beta }_{123}$	2	2.442	2.115	2.094	1.565	0.823	0.638
13	${\Delta }_{11}$	0.1	0.636	0.44	0.155	5.038	1.941	1.354
14	${\Delta }_{12}$	0.15	0.677	0.414	0.269	5.101	1.98	1.387
15	${\Delta }_{21}$	0.25	0.587	0.461	0.25	5.047	1.952	1.354
16	${\Delta }_{22}$	0.3	0.671	0.455	0.351	5.145	1.978	1.383
17	${\gamma }_{11}$	−1.5	−1.926	−1.585	−1.662	5.088	1.992	1.396
18	${\gamma }_{12}$	2	0.475	0.226	0.929	26.556	9.533	6.397
19	${\gamma }_{22}$	3	7.729	6.204	5.216	40.777	13.805	8.713
20	${\gamma }_{13}$	2	2.33	2.134	2.273	5.22	2.008	1.447

. The results show that the estimated parameters approach their true coefficient on average and the corresponding standard errors also decrease as the sample size increases.

Furthermore, the suitability of the estimated parameter with the true coefficients for each sample size can be presented in Figure 1 below:

Based on Figure 1, the results show that the greater of samples size, the closer the estimated parameter values get to the true coefficient values. This can be seen from the movement of the estimated parameter coefficients for n = 300, which are closest to the true coefficients, and for n = 100, the coefficients which are farthest for each parameter. In addition, the standard error range values from parameter estimation was also investigated for different samples in the form of boxplots. When the sample size increases, the standard error values of the model parameter estimation becomes less, as presented in Figure 2 as follows:

6. Conclusions

We developed ordinal logistic regression with two correlated response variables in which the relationship between the continuous predictor variable and its logit is modeled as a polynomial form, called the BPOLR model. The proposed BPOLR model contributes to the development of statistical science that can be used as an alternative solution in statistical modeling related to the bivariate ordinal logistic regression, where the relationship between continuous predictor variables and their logit links is assumed to be not always linear, so that it is more flexible. Therefore, the advantage of the BPOLR model can accommodate linear, quadratic, cubic or other relationships between continuous predictor variables and their logit links in a model. We used the MLE method to obtain parameter estimators of the BPOLR model. The first partial derivatives are not closed-form, thus, a numerical optimization such as the method of BHHH is needed to obtain the maximum likelihood estimators. We used the MLRT method to obtain test statistics of parameters. The simultaneous test was used to make the simultaneous of the parameters significantly. The statistical test distribution has a Chi-square limit that, asymptotically with the degree of freedom, is the difference between the parameter effective number in the reduced and full models. The optimal degree value for each predictor variable is obtained when the value of AICc, BIC, and deviance are minimum. Based on the simulation study results, the estimated parameters approach their true coefficient on average and the corresponding standard errors also decrease as the sample size increases.