Dichotomous Proportional Hazard Regression Model: A Case Study on Students’ Dropout

Abstract

In problems involving binary classification, researchers often encounter data suitable for modeling dichotomous responses. These scenarios include medical diagnostics, where outcomes are classified as “disease” or “no disease”, and credit scoring in finance, determining whether a loan applicant is “high risk” or “low risk”. Dichotomous response models are also useful in many other areas for estimating binary responses. The logistic regression model is one option for modeling dichotomous responses; however, other statistical models may be required to improve the quality of fits. In this paper, a new regression model is proposed for cases where the response variable is dichotomous. This novel, non-linear model is derived from the cumulative distribution function of the proportional hazard distribution, and is suitable for modeling binary responses. Statistical inference is performed using a classical approach with the maximum likelihood method for the proposed model. Additionally, it is demonstrated that the introduced model has a non-singular information matrix. The results of a simulation study, along with an application to student dropout data, show the great potential of the proposed model in practical and everyday situations.

1. Introduction

In recent statistical literature, new probability distributions have been introduced as extensions of other known distributions. This methodology serves as a foundation for generating new families of distributions applicable across various fields. Among other authors, this approach was utilized by Eugene et al. [1] to propose the Beta-G class of distributions. Subsequently, Silva et al. [2] introduced the modified Weibull beta distribution families and the Weibull beta geometry, as noted by Cordeiro et al. [3]. Moreover, building upon this methodology, Cordeiro and de Castro [4] defined the Kumaraswamy-G class of distributions, followed by the suggestion of the Kumaraswamy modified Weibull by Cordeiro et al. [5]. Similarly, Zografos and Balakrishnan [6] and Ristić and Balakrishnan [7] presented a new family of distributions generated by gamma random variables, leading to the development of the Gamma-Generated-Logistic distributions by Castellares et al. [8] and the Gamma-Birnbaum–Saunders distributions by Cordeiro et al. [9].

On the other hand, Martínez-Flórez et al. [10] examined the exponentiated-skew-normal distribution. Similarly, Martínez-Flórez et al. [11] proposed the family of proportional hazard distributions based on the distribution of the minimum in the sample. All of these new families of distributions have proven useful in analyzing responses of interest by adjusting both linear and nonlinear regression models. For instance, the regression model with skew-normal distributed errors Azzalini [12] has been widely utilized. Furthermore, extensions of regression models have been recommended, assuming errors follow the exposed skew-normal distributions (Martínez-Flórez et al. [10]; Martínez-Flórez et al. [13]). Moreover, these distribution families have been extended to encompass the case of the Birnbaum distribution [14] and the Birnbaum–Saunders log-linear regression model proposed by Rieck and Nedelman [15], showcasing the extensive range of symmetric and asymmetric families available in the literature.

All the works previously presented, and many others that have not been mentioned here, are appropriate in cases where the response variable has its support in the set of real numbers or has positive support, while very few works focus on the problem of dealing with dichotomous data. In this particular case, the issue is addressed based on non-linear functions or link functions such as the logistic regression model, known in the statistical literature as the logit model, or the non-linear alternative based on the cumulative distribution function (CDF) of the normal density, called the probit model. Thus, the limited existence of proposals in the statistical literature for the analysis of dichotomous or polytomous responses through link functions used in other types of models becomes evident.

In practice, the regression model with a dichotomous response (logistic model) has been widely used in several areas of knowledge. In the educational area, for example, it can be used to predict the probability of a student dropping out based on their academic performance, age of entry, the educational level of their parents, number of siblings, etc. In the health sector, certain patient characteristics and the application of specific treatments can be analyzed using the model to understand the connection between the patient and the implemented treatment, including the probability or odds of survival. Similarly, in finance, based on characteristics such as sex, age, race, income, and educational level, the behavior of investors can be predicted. These models are also utilized for classifying individuals into certain groups according to the predicted probability of a specific event occurring.

In this article, a new regression model is proposed to address research with dichotomous response variables. This novel model can be applied to various fields, including medicine, finance, education, and the social sciences. Our proposal is grounded in the family of proportional hazard distributions, specifically utilizing an extension of the logistic distribution within this family.

The remainder of this work is organized as follows: Section 2 provides a brief description of the logistic distribution and its associated regression model. Section 3 describes the proportional hazard and proportional hazard logistic distributions, along with some of their most important properties. In Section 4, the proportional hazard logistic regression model is introduced. Additionally, the statistical inference process is performed using a classical approach, presenting the score function and the elements of the observed information matrix. Section 5 presents an application of the introduced model to student dropout data. Finally, the conclusions of the paper are presented in Section 6.

2. Logistic Distribution

A continuous random variable with a logistic distribution has a probability density function (PDF) given by

$$f_{\mathrm{L}}\left(z\right)=\frac{exp(-z)}{{(1+exp(-z))}^2}=\frac{1}{4}{\mathrm{sech}}^2\left(\frac{z}{2}\right),z\in \mathbb{R}.$$

(1)

where $\mathrm{sech}$ denotes the hyperbolic secant function. The shape of the logistic distribution is similar to the shape of the normal density, with heavier tails and greater kurtosis than the normal distribution.

The cumulative distribution function (CDF) of a random variable with a logistic distribution is given by

$${\mathcal{F}}_{\mathrm{L}}\left(z\right)=\frac{exp\left(z\right)}{1+exp\left(z\right)}=\frac{1}{2}+\frac{1}{2}tanh\left(\frac{z}{2}\right),$$

while its survival and hazard functions can be written as

$${\mathcal{S}}_{\mathrm{L}}\left(z\right)=\frac{1}{1+exp\left(z\right)}=\frac{1}{2}-\frac{1}{2}tanh\left(\frac{z}{2}\right)\mathrm{and}$$

$$h_{\mathrm{L}}\left(z\right)=\frac{exp\left(z\right)}{1+exp\left(z\right)},$$

where tanh denotes the hyperbolic tangent function.

The extension of the logistic distribution to the location-scale case is achieved by using the transformation $Y=\mu +\sigma Z$ with $\mu \in \mathbb{R}$ and $\sigma >0$. This is denoted by $Y\sim \mathrm{L}(\mu ,\sigma )$, where $\mu$ represents the location parameter and $\sigma$ the scale. Since this distribution is symmetric, then $\mathbb{E}\left(Y\right)=\mu$, $\mathrm{Var}\left(Y\right)=\frac{{\pi }^2}{3}{\sigma }^2$, the asymmetry coefficient is zero, and its excess kurtosis is equal to $\frac{6}{5}$. Finally, the p-th percentile, for $0<p<1$, of this distribution is given by $y_p=\mu +\sigma log(p/(1-p))$.

Associated with the logistic distribution is the logistic regression model, which is used to explain the probability of success of a random variable with a binomial distribution when there is a set of covariates that explain this probability (see Agresti [16]). In essence, the logistic regression model is given by

$$p_i=\mathrm{Pr}(Y_i=1\mid x_1,x_2,\dots ,x_p)=\frac{exp\left({\mathit{x}}_i^{\top }\mathit{\beta }\right)}{1+exp\left({\mathit{x}}_i^{\top }\mathit{\beta }\right)}$$

where $\mathit{x}={(1,x_1,\dots ,x_p)}^{\top }$ represents a vector of covariates, $\mathit{\beta }={({\beta }_0,{\beta }_1,\dots ,{\beta }_p)}^{\top }$ is the vector of model coefficients (unknown values that must be estimated), and $Y_i$ is a Bernoulli random variable with parameter $p_i$.

3. Hazard Proportional Distribution

In recent decades, families of asymmetric distributions have been introduced for fitting data with tails heavier or lighter than the normal distribution. As is well known, in the presence of high degrees of skewness and/or kurtosis, inferential processes based on the assumption of normality are inadequate. Similarly, while the elliptical family may provide a solution for distributions with heavy tails, it fails to address the issue of asymmetry in the data under study.

The skew-normal (SN) distribution, introduced by Azzalini [12], is defined by the PDF given as

$$\phi (z;\lambda )=2\varphi \left(z\right)\Phi \left(\lambda z\right),z\in \mathbb{R},$$

(2)

where $\varphi$ and $\Phi$ represent the PDF and CDF of the standard normal distribution, respectively, and $\lambda$ is a skewness parameter. The distribution is denoted by $Z\sim SN\left(\lambda \right)$. In addition to the work of Azzalini [12], the SN distribution described in (2) has been extensively studied by Henze [17], Pewsey [18], Chiogna [19], and Gómez et al. [20], among others.

Building on the work of Lehmann [21], Martínez-Flórez et al. [11] investigated another family of asymmetric univariate distributions called the proportional hazard. The PDF of this distribution is given by

$${\phi }_{\mathrm{F}}(z;\alpha )=\alpha f\left(z\right){\{1-F\left(z\right)\}}^{\alpha -1},z\in \mathbb{R},$$

(3)

where $\alpha$ is a positive real number, and F is a continuous CDF with continuous PDF f. This distribution is denoted by $\mathrm{PHF}\left(\alpha \right)$. The hazard function associated with the density ${\phi }_F$ is

$$h_{{\phi }_F}(X,\alpha )=\alpha h_f\left(x\right),$$

where $h_f=f/(1-F)$ represents the hazard function related to the density f. When $F=\Phi (·)$ and $f=\varphi (·)$, the distribution is called proportional hazard normal, denoted by $\mathrm{PHN}\left(\alpha \right)$. The PDF is given by

$${\phi }_{\Phi }(z;\alpha )=\alpha \varphi \left(z\right){\left\{\mathcal{S}\left(z\right)\right\}}^{\alpha -1},z\in \mathbb{R},$$

(4)

where $\mathcal{S}\left(z\right)$ is the survival function associated with the PDF $\varphi (·)$. This model serves as an alternative to accommodate data with asymmetry and kurtosis that fall outside the ranges allowed by the normal distribution.

The CDF of the $\mathrm{PHN}\left(\alpha \right)$ distribution is given by:

$$F_{\Phi }(z;\alpha )=1-{\left\{\mathcal{S}\left(z\right)\right\}}^{\alpha },z\in \mathbb{R}.$$

(5)

By varying the $\alpha$ parameter, Martínez-Flórez et al. [11] found that the range of asymmetry and kurtosis coefficients, $\sqrt{{\beta }_1}$ and ${\beta }_2$, respectively, of the variable Z∼PHN($\alpha$) falls within the intervals $(-1.1578,0.9918)$ and $(1.1513,4.3023)$. These ranges exhibit better skewness and kurtosis properties than those of the SN distribution. Additionally, Martínez-Flórez et al. [11] demonstrated that the information matrix of the PHN distribution in the location-scale case, denoted as $\mathrm{PHN}(\mu ,\sigma ,\alpha )$, is nonsingular. A particular case of the proportional hazard family is discussed below.

Proportional Hazard Logistic Distribution

The proportional hazard logistic (PHL) distribution, denoted by $\mathrm{PHL}\left(\alpha \right)$, is defined by the PDF given as

$$\begin{array}{ccc}\hfill {\phi }_{\mathrm{HL}}\left(x;\alpha \right)& =& \alpha \frac{exp\left(x\right)}{{\left(1+exp\left(x\right)\right)}^{\alpha +1}}\hfill \\ & =& \frac{\alpha }{4}{\mathrm{sech}}^2\left(\frac{x}{2}\right){\left[\frac{1}{2}-\frac{1}{2}tanh\left(\frac{x}{2}\right)\right]}^{\alpha -1}\hfill \end{array}$$

(6)

Its respective CDF is given by

$$\begin{array}{ccc}\hfill {\mathcal{F}}_{\mathrm{HL}}\left(x;\alpha \right)& =& 1-\frac{1}{{\left(1+exp\left(x\right)\right)}^{\alpha }}\hfill \\ & =& 1-{\left[\frac{1}{2}-\frac{1}{2}tanh\left(\frac{x}{2}\right)\right]}^{\alpha }\hfill \end{array}$$

(7)

while the survival and hazard functions can be expressed as

$${\mathcal{S}}_{\mathrm{HL}}\left(x;\alpha \right)=\frac{1}{{\left(1+exp\left(x\right)\right)}^{\alpha }}={\left[\frac{1}{2}-\frac{1}{2}tanh\left(\frac{x}{2}\right)\right]}^{\alpha }$$

(8)

and

$$\begin{array}{ccc}\hfill h_{\mathrm{HL}}\left(x;\alpha \right)& =& \alpha \frac{exp\left(x\right)}{1+exp\left(x\right)}\hfill \\ & =& \frac{\alpha }{4}\frac{{\mathrm{sech}}^2\left(\frac{x}{2}\right)}{\frac{1}{2}-\frac{1}{2}tanh\left(\frac{x}{2}\right)}=\alpha h_L\left(x\right).\hfill \end{array}$$

(9)

respectively, where $h_L\left(x\right)$ is the hazard function of the logistic distribution.

Figure 1 illustrates the behavior of the CDF and the survival function for different values of the parameter $\alpha$. It is noteworthy that, for $\alpha =1$, the CDF corresponds to that of the logistic distribution. Moreover, the hazard function of the PHL is a multiple of the hazard function of the logistic distribution. Additionally, the adjustment of the CDF of the PHL distribution is more flexible than that of the logistic distribution. Similarly, it is observed that for $\alpha =0.75$, the survival function converges more slowly (indicating a higher probability of survival) towards zero compared to the survival function of the logistic distribution, whereas for values greater than zero, the convergence to zero is faster (indicating a lower probability of survival).

The r-th moment of the random variable $Y\sim \mathrm{PHL}\left(\alpha \right)$ is given by:

$$\mathbb{E}\left(Y^r\right)={\int }_1^{\infty }\frac{{log}^r(u-1)}{{(u-1)}^2}{(1-u^{-1})}^{\alpha +1}du.$$

(10)

From Expression (10), the moments of orders 1, 2, 3, and 4 of the PHL distribution can be derived, facilitating the numerical calculation of its mean, variance, skewness, and kurtosis coefficients.

4. Proportional Hazard Logistic Regression Model

Assuming the regression model:

$$Y_i={\mathit{X}}_i^{\top }\mathit{\beta }+{\epsilon }_i={\mu }_i+{\epsilon }_i,i=1,2,\dots ,n$$

(11)

where $\mathit{X}={(1,x_1,\dots ,x_p)}^{\top }$ represents a set of covariates, $\mathit{\beta }={({\beta }_0,{\beta }_1,\dots ,{\beta }_p)}^{\top }$ denotes a set of unknown coefficients, and ${\epsilon }_i$∼$\mathrm{PHL}(0,\sigma ,\alpha )$. It then follows that Y_i∼PHL$({\mu }_i,\sigma ,\alpha )$, for $i=1,2,\dots ,n$.

However, when Y is a dichotomous random variable with values zero and one, the model errors are not independent and do not satisfy the assumption of homoscedasticity. Additionally, it cannot be ensured that

$$\mathbb{E}(Y_i\mid x_1,\dots ,x_p)=\mathrm{Pr}(Y_i=1\mid x_1,\dots ,x_p)$$

is bounded by 0 and 1.

For this reason, it is necessary to determine a distribution function $G(·)$ such that

$$\mathrm{Pr}(Y_i=1\mid x_1,\dots ,x_p)=p_i=G(Y_i=1\mid x_1,\dots ,x_p).$$

The function $G(·\mid x_1,\dots ,x_p)$ is known as a link function, and since it must ensure that the prediction lies between 0 and 1, it is commonly chosen as the distribution function of certain random variables studied in classical probability theory literature.

The link functions $G(·\mid x_1,\dots ,x_p)$ typically utilized in practice are the CDF of the logistic distribution, resulting in the logit model, and the CDF of the normal distribution, resulting in the probit model. Due to their mathematical and computational complexity, the logit model is generally preferred over the probit model in practical applications. A notable commonality between these two models is their symmetric CDF, which can be a limitation in scenarios where the probability of success for response variable Y exhibits asymmetric behavior. Moreover, both distributions have limitations in accurately modeling certain probabilities in their tails. As illustrated in Figure 1, the CDF of the proportional hazard logistic distribution displays asymmetric behavior. Additionally, the inclusion of the parameter $\alpha$ allows for modeling the probabilities in its tails. This parameter enhances the flexibility of the probability of success compared to the logit and probit functions, suggesting the potential for more precise adjustment of the probability of success for the variables under study.

Referring to $G(·\mid x_1,\dots ,x_p)$ as the CDF of the PHL, it follows that

$$\begin{array}{cc}\hfill \mathrm{Pr}(Y_i=1\mid x_1,\dots ,x_p)& =p_i=G(Y_i=1\mid x_1,\dots ,x_p)\hfill \\ \hfill & =1-\frac{1}{{\left(1+exp\left({\mathit{x}}_i^{\top }\mathit{\beta }\right)\right)}^{\alpha }}\hfill \\ \hfill & =1-{\left[\frac{1}{2}-\frac{1}{2}tanh\left(\frac{{\mathit{x}}_i^{\top }\mathit{\beta }}{2}\right)\right]}^{\alpha }.\hfill \end{array}$$

From this, it is obtained that

$$\begin{array}{cc}\hfill \mathrm{Pr}(Y_i=0\mid x_1,\dots ,x_p)& =1-\mathrm{Pr}(Y_i=1\mid x_1,\dots ,x_p)\hfill \\ \hfill & =\frac{1}{{\left(1+exp\left({\mathit{x}}_i^{\top }\mathit{\beta }\right)\right)}^{\alpha }}\hfill \\ \hfill & ={\left[\frac{1}{2}-\frac{1}{2}tanh\left(\frac{{\mathit{x}}_i^{\top }\mathit{\beta }}{2}\right)\right]}^{\alpha }.\hfill \end{array}$$

For $p_i=\mathrm{Pr}(Y_i=1\mid x_1,\dots ,x_p)$, it follows that

$$\begin{array}{c}\hfill log\left(\frac{1-{(1-p_i)}^{1/\alpha }}{{(1-p_i)}^{1/\alpha }}\right)={\mathit{x}}_i^{\top }\mathit{\beta },i=1,2,\dots ,n,\end{array}$$

(12)

which will be referred to as the logit complement $\alpha -$root transformation.

4.1. Properties of the PHL Regression Model

Given the structure of the probability function included in this new model, some statistics of interest are calculated for the interpretation of the parameters. Then, the odds $\mathrm{odds}(x_1,x_2,x_3,\dots ,x_p)=\mathrm{odds}\left(\mathit{x}\right)$ are given by

$$\begin{array}{ccc}\hfill \mathrm{odds}\left({\mathit{x}}_i\right)& =& \frac{\mathrm{Pr}(Y_i=1\mid x_1,x_2,\dots ,x_p)}{1-\mathrm{Pr}(Y_i=1\mid x_1,x_2,\dots ,x_p)}\hfill \\ & =& {\left(1+exp\left({\mathit{x}}_i^{\top }\mathit{\beta }\right)\right)}^{\alpha }-1.\hfill \end{array}$$

Thus, the relative risk $\left(RR\right)$ or odds ratio, to compare the profile of individuals i and k, is given by

$$RR(i,k)=\frac{\mathrm{odds}\left({\mathit{x}}_i\right)}{\mathrm{odds}\left({\mathit{x}}_k\right)}=\frac{{\left(1+exp\left({\mathit{x}}_i^{\top }\mathit{\beta }\right)\right)}^{\alpha }-1}{{\left(1+exp\left({\mathit{x}}_k^{\top }\mathit{\beta }\right)\right)}^{\alpha }-1}.$$

This expression is used when there are profiles of different individuals, or when the profiles only differ in the jth variable. Thus, to estimate the relative risk in the ith individual when the jth variable is increased by one unit, denoted as $x_j+1$, while keeping the value of the rest of the variables constant, we have the expression

$$\frac{\mathrm{odds}(x_1,\dots ,x_{j-1},x_j+1,x_{j+1},\dots ,x_p)}{\mathrm{odds}(x_1,\dots ,x_{j-1},x_j,x_{j+1},\dots ,x_p)}=$$

$$\frac{{\left(1+exp\left({\beta }_j\right)exp\left({\mathit{x}}_i^{\top }\mathit{\beta }\right)\right)}^{\alpha }-1}{{\left(1+exp\left({\mathit{x}}_i^{\top }\mathit{\beta }\right)\right)}^{\alpha }-1}$$

This represents the odds or the number of times the risk of the event occurring increases (or decreases) when the variable $x_j$ increases by one unit.

4.2. Maximum Likelihood Estimation

Given a random sample $y_1,y_2,\dots ,y_n$ of a random variable Y with distribution $Y_i$∼$\mathrm{Bin}(n,p_i)$, and considering a set of covariates $x_1,x_2,\dots ,x_p$, the likelihood function is expressed as

$${\mathcal{L}}_{\mathrm{PHL}}(\mathit{\beta },\alpha \mid \mathit{X},\mathit{Y})=\prod _{i=1}^n{p}_i^{y_i}{(1-p_i)}^{1-y_i}.$$

Then, the log-likelihood function is given by

$$\begin{array}{ccc}\hfill {\ell }_{\mathrm{PHL}}(\mathit{\beta },\alpha |\mathit{X},\mathit{Y})& =& \sum _{i=1}^n{y}_{i}log\left(p_i\right)+(1-y_i)log(1-p_i)\hfill \\ & =& \sum _{i=1}^n{y}_{i}log\left({\left(1+exp\left({\mathit{x}}_i^{\top }\mathit{\beta }\right)\right)}^{\alpha }-1\right)\hfill \\ & & -\alpha \sum _{i=1}^{n}log\left(1+exp\left({\mathit{x}}_i^{\top }\mathit{\beta }\right)\right)\hfill \end{array}$$

(13)

The score function, $U(\mathit{\beta },\alpha )=\left(U\right(\mathit{\beta }),U(\alpha \left)\right)$ with $U\left(\mathit{\beta }\right)=(U\left({\beta }_0\right),U\left({\beta }_1\right),U\left({\beta }_2\right),\dots ,$ $U\left({\beta }_p\right))$, which is calculated as the first derivative of the log-likelihood function concerning the parameters, is given by

$$\begin{array}{ccc}\hfill U\left({\beta }_j\right)& =& \alpha \sum _{i=1}^n{x}_{ij}{y}_{i}exp\left({\mathit{x}}_i^{\top }\mathit{\beta }\right)\frac{{\left(1+exp\left({\mathit{x}}_i^{\top }\mathit{\beta }\right)\right)}^{\alpha -1}}{{\left(1+exp\left({\mathit{x}}_i^{\top }\mathit{\beta }\right)\right)}^{\alpha }-1}\hfill \\ & +& \alpha \sum _{i=1}^n{x}_{ij}\frac{exp\left({\mathit{x}}_i^{\top }\mathit{\beta }\right)}{1+exp\left({\mathit{x}}_i^{\top }\mathit{\beta }\right)}\hfill \end{array}$$

(14)

for $j=0,1,2,\dots ,p$, and

$$\begin{array}{ccc}\hfill U\left(\alpha \right)& =& \sum _{i=1}^n{y}_i\frac{{\left(1+exp\left({\mathit{x}}_i^{\top }\mathit{\beta }\right)\right)}^{\alpha }log\left(1+exp\left({\mathit{x}}_i^{\top }\mathit{\beta }\right)\right)}{{\left(1+exp\left({\mathit{x}}_i^{\top }\mathit{\beta }\right)\right)}^{\alpha }-1}\hfill \\ & -& \sum _{i=1}^{n}log\left(1+exp\left({\mathit{x}}_i^{\top }\mathit{\beta }\right)\right)\hfill \end{array}$$

(15)

The elements of the observed information matrix, $\kappa \left(\mathit{\theta }\right)$, defined as minus the Hessian matrix (matrix of second derivatives concerning the parameters), are given by:

$$\begin{array}{ccc}\hfill {\kappa }_{{\beta }_j{\beta }_k}& =& \alpha \sum _{i=1}^n{x}_{ij}{x}_{ik}\frac{exp\left({\mathit{x}}_i^{\top }\mathit{\beta }\right)}{{\left(1+exp\left({\mathit{x}}_i^{\top }\mathit{\beta }\right)\right)}^2}\left[1+\frac{y_i}{p_i^2}\right.\hfill \\ & & \left.\left(-p_i\left(1+exp\left({\mathit{x}}_i^{\top }\mathit{\beta }\right)\right)+exp\left({\mathit{x}}_i^{\top }\mathit{\beta }\right)(p_i+\alpha (1-p_i))\right)\right]\hfill \\ \hfill {\kappa }_{{\beta }_j\alpha }& =& \sum _{i=1}^n{x}_{ij}\frac{exp\left({\mathit{x}}_i^{\top }\mathit{\beta }\right)}{1+exp\left({\mathit{x}}_i^{\top }\mathit{\beta }\right)}\left[1-\frac{y_i}{p_i^2}\left(\frac{p_i}{1+exp\left({\mathit{x}}_i^{\top }\mathit{\beta }\right)}-\right.\right.\hfill \\ & & \left.\left.\alpha log(1+exp\left({\mathit{x}}_i^{\top }\mathit{\beta }\right))(1-p_i)\right)\right]\hfill \\ \hfill {\kappa }_{\alpha \alpha }& =& \sum _{i=1}^n\frac{1-p_i}{p_i^2}{log}^2\left(1+exp\left({\mathit{x}}_i^{\top }\mathit{\beta }\right)\right).\hfill \end{array}$$

The elements of the information matrix, which are obtained from the expected value of the elements of the observed information matrix, $I\left(\mathit{\theta }\right)=\mathbb{E}\left(\kappa \left(\mathit{\theta }\right)\right)$, are given by

$$\begin{array}{ccc}\hfill i_{{\beta }_j{\beta }_k}& =& \alpha \sum _{i=1}^n{x}_{ij}{x}_{ik}\frac{1-p_i}{p_i}{\left(\frac{exp\left({\mathit{x}}_i^{\top }\mathit{\beta }\right)}{1+exp\left({\mathit{x}}_i^{\top }\mathit{\beta }\right)}\right)}^2,\hfill \\ \hfill i_{{\beta }_j\alpha }& =& \sum _{i=1}^n{x}_{ij}{\left(\frac{exp\left({\mathit{x}}_i^{\top }\mathit{\beta }\right)}{1+exp\left({\mathit{x}}_i^{\top }\mathit{\beta }\right)}\right)}^2-\hfill \\ & & \alpha \sum _{i=1}^n{x}_{ij}\frac{1-p_i}{p_i}\frac{exp\left({\mathit{x}}_i^{\top }\mathit{\beta }\right)}{1+exp\left({\mathit{x}}_i^{\top }\mathit{\beta }\right)}\hfill \\ & & log\left(1+exp\left({\mathit{x}}_i^{\top }\mathit{\beta }\right)\right)\hfill \\ \hfill i_{\alpha \alpha }& =& \sum _{i=1}^n\frac{1-p_i}{p_i^2}{log}^2\left(1+exp\left({\mathit{x}}_i^{\top }\mathit{\beta }\right)\right).\hfill \end{array}$$

When $\alpha =1$, we obtain $p_i=\frac{exp\left({\mathit{x}}_i^{\top }\mathit{\beta }\right)}{1+exp\left({\mathit{x}}_i^{\top }\mathit{\beta }\right)}$, and the information matrix can be written as

$$I_F\left(\theta \right)=\left(\begin{array}{cc}{\mathit{X}}^{\top }\mathit{W}\mathit{X}& {\mathit{X}}^{\top }\mathit{W}\\ {\mathit{M}}^{\top }\mathit{X}& \frac{1-p_i}{p_i^2}{log}^2(1+exp\left({\mathit{x}}_i^{\top }\mathit{\beta }\right))\end{array}\right),$$

(16)

where $\mathit{W}$ is the diagonal matrix $\mathit{W}=\mathrm{diag}\left(p_i(1-p_i)\right),i=1,2,\dots ,n$, and $\mathit{M}$ is a vector with elements $m_i=p_i\left(p_i-\frac{1-p_i}{p_i}log(1+exp\left({\mathit{x}}_i^{\top }\mathit{\beta }\right))\right)$.

Letting $d=\frac{1-p_i}{p_i^2}{log}^2(1+exp\left({\mathit{x}}_i^{\top }\mathit{\beta }\right))$, we obtain that the determinant of the information matrix is given by:

$$\left|I\left(\mathit{\theta }\right)\right|=d^{-p}\left|{\mathit{X}}^{\top }\left(\mathit{W}-k^{-1}\mathit{M}{\mathit{M}}^{\top }\right)\mathit{X}\right|\ne 0.$$

Thus, the information matrix is non-singular, which guarantees the existence of the variance–covariance matrix of the vector of maximum likelihood estimators (MLE) $\widehat{\mathit{\theta }}$. It can also be concluded that the variance–covariance matrix of the MLE can be written as:

$$\mathbf{\Sigma }=I^{-1}\left(\widehat{\mathit{\theta }}\right).$$

Therefore, for large sample sizes, we have

$$\widehat{\mathit{\theta }}\stackrel{d}{\to }{N}_{p+2}(\mathit{\theta },\Sigma ),$$

meaning that the distribution of the vector of estimators is consistent and asymptotically normal, with a covariance matrix equal to the inverse of the Fisher information matrix.

Confidence intervals for coefficients ${\theta }_r$ of level $100(1-\psi )\%$ can be obtained from the expression ${\widehat{\theta }}_r\mp z_{1-\psi /2}\sqrt{\widehat{\sigma }\left({\widehat{\theta }}_r\right)}$. Additionally, the adequacy of the proportional hazard logistic regression (PHLR) model can be evaluated through hypothesis testing:

$$H_0:{\beta }_1={\beta }_2=\dots ={\beta }_p=0\mathrm{vs}{H}_1:{\beta }_j\ne 0,$$

for at least one $j=1,\dots ,p$.

We can use the deviance function given by

$$G_p=-2(\ell ({\beta }_0,\alpha )-\ell (\widehat{\mathit{\beta }},\alpha )),$$

with distribution $G_p\sim {\chi }_p^2$. Similarly, two models can be compared: one complete with r variables $\left({\mathit{\beta }}_r\right)$, and another with q$(q<r)$ variables $\left({\mathit{\beta }}_q\right)$ through the test statistic

$$G_{r-q}=-2(\ell ({\widehat{\mathit{\beta }}}_q,\alpha )-\ell ({\widehat{\mathit{\beta }}}_r,\alpha )),$$

for which we have $G_{r-q}\sim {\chi }_{r-q}^2$. This same statistic is useful to test the significance of the remaining $r-q$ variables that were not included in the model with q variables.

One of the strategies to validate the good fit of the logistic regression model is to analyze the proportion of correct classification that the fitted model achieves. Letting $G_1$ be the group of observations with $Y_i=1$, and $G_2$ be the group of observations with $Y_i=0$ then, using Bayes’ Theorem, the probability of classifying an individual into group $G_1$ given the information of the explanatory variables $x_1,x_2,\dots ,x_p$ is given by

$$\mathrm{Pr}(G_1\mid \mathit{x})=\frac{p_1\times \mathrm{Pr}(\mathit{x}\mid G_1)}{p_1\times \mathrm{Pr}(\mathit{x}\mid G_1)+p_2\times \mathrm{Pr}(\mathit{x}\mid G_2)}.$$

Thus, when performing the calculations for our model, we have

$$\mathrm{Pr}(G_1\mid \mathbf{x})=1-\frac{p_2}{(p_2-p_1)+p_1{(1+exp\left({\mathit{x}}_i^{\top }\mathit{\beta }\right))}^{\alpha }}.$$

Similarly, $\mathrm{Pr}(G_2\mid \mathit{x})$ is defined. In this case, the decision will be to classify the ith individual into $G_1$ if $\mathrm{Pr}(G_1\mid \mathit{x})>\mathrm{Pr}(G_2\mid \mathit{x})$; that is, if $\mathrm{Pr}(G_1\mid \mathit{x})>0.5$.

To evaluate the predictive capacity of the proportional hazard logistic regression model, the overall accuracy of the model can be calculated, which is defined as the proportion of individuals that are correctly classified, as well as the sensitivity or true positive rate of the model (TPR), defined as the number of correctly classified individuals from group $G_1$ divided by the total number of correctly classified ones (from $G_1$ and $G_2$). Similarly, the false negative rate (FNR) of the model is defined as $(1-TPR)$, among others.

5. Case Study: Students’ Dropout Data

The data for this application consist of a sample of 413 students from the Department of Mathematics and Statistics of the University of Córdoba, which were obtained from the SPADIES System of the Ministry of National Education of Colombia (MNE). The response variable in this application takes the values $Y=1$ (if program dropout) or $Y=0$ (if non-dropout). The explanatory variables considered are $x_1=$ (cumulative general average, CGA), $x_2=$ character of the school (CS) of the student where they studied, taking values $=1$ (if the student comes from an official school), and $=0$ (if not), and $x_3=$ the number of periods enrolled (NPE). The logistic regression (LR) and proportional hazard logistic regression (PHLR) models were fitted. The results of the fitted models, obtained using the R Development Core Team [22] package, are given in

Table 1. Parameter estimation of LR and PHLR models (standard errors of the estimates are given in parentheses).

Model	${\widehat{\beta }}_0$	${\widehat{\beta }}_1$	${\widehat{\beta }}_2$	${\widehat{\beta }}_3$	$\widehat{\alpha }$	AIC	CAIC	BIC
LR	7.290	−1.576	0.442	−0.270		310.4	330.5	326.5
se	(1.220)	(0.384)	(0.423)	(0.035)
PHLR	1.374	−0.860	0.355	−0.201	11.93	295.9	325.0	316.0
se	(1.206)	(0.340)	(0.293)	(0.025)	(5.301)

.

The results of the model fit indicate that the variables CGA and NPE are significant, whereas the variable CS does not significantly explain the probability of university student dropout.

To compare the fitted models, we employ the Akaike Information Criterion (AIC) Akaike [23], corrected AIC (CAIC), and the Bayesian Information Criterion (BIC) by Hastie and Tibshirani [24], given by

$$AIC=-2\times \widehat{\ell }(·)+2p,$$

$$CAIC=-2\times \widehat{\ell }(·)+2p\left(1+\frac{n+2}{n-p-2}\right)$$

and

$$BIC=-2\times \widehat{\ell }(·)+plog\left(n\right)$$

where p is the number of parameters in the model and n is the sample size. The results favor the PHLR model based on AIC, CAIC, and BIC values.

To compare the proportional hazard logistic regression (PHLR) model with the logistic regression model, we conduct the hypothesis test

$$H_0:\alpha =1\mathrm{vs}{H}_1:\alpha \ne 1,$$

using the likelihood ratio statistic

$${\Lambda }_1=\frac{{\mathcal{L}}_{\mathrm{L}}\left(\mathit{\theta }\right)}{{\mathcal{L}}_{\mathrm{PHL}}\left({\mathit{\theta }}^{\ast }\right)},$$

where ${\mathcal{L}}_{\mathrm{L}}(·)$ and ${\mathcal{L}}_{\mathrm{PHL}}(·)$ represent the likelihood functions of the logistic and PHL models, respectively. Upon numerical evaluation, we obtain

$$-2log(\Lambda )=-2(-151.2+142.955)=16.49,$$

which exceeds the value of ${\chi }_{1,95\%}^2=3.84$. The PHL model exhibits the best fit compared to the logistic model.

Carrying out the hypothesis test of the significance of the explanatory variables

$$H_0:{\beta }_1={\beta }_2={\beta }_3=0\mathrm{vs}{H}_1:{\beta }_j\ne 0,$$

for at least one $j=1,2,3$, we have

$$\begin{array}{ccc}\hfill G_3^2& =& -2(-262.4158+142.9585)\hfill \\ & =& 238.9147>{\chi }_{0.05,3}^2=7.8147,\hfill \end{array}$$

therefore, the null hypothesis is rejected. Similarly, for the hypothesis test

$$H_0:{\beta }_2=0\mathrm{vs}{H}_1:{\beta }_2\ne 0,$$

it follows that

$$\begin{array}{cc}\hfill G_3^2& =-2(\ell ({\beta }_0,{\beta }_1,{\beta }_3,\alpha )-\ell ({\beta }_0,{\beta }_1,{\beta }_2,{\beta }_3,\alpha ))\hfill \\ \hfill & =-2(-143.96+142.9585)=2.003<{\chi }_{0.05,1}^2=3.84,\hfill \end{array}$$

Therefore, the null hypothesis is not rejected, meaning the variable character of the school is not significant in the model. However, academic differences are observed in the classroom between students who come from official schools and those from private schools, with the latter demonstrating better preparation.

Note that in the proportional hazard logistic regression model, the case $\alpha =1$ corresponds to the logistic distribution. However, the hypothesis test $H_0:\alpha =1$ vs. $H_1:\alpha \ne 1$, which is performed using the likelihood ratio statistic, is rejected. This means that the parameter $\alpha$ is significantly different from one, and must be considered to explain the behavior of the data. Moreover, the AIC, CAIC, and BIC criteria are favorable to the logistic proportional hazard model when compared with the usual logistic regression model. All of the above allows us to conclude that the proportional hazard logistic regression model fits better.

So, the fitted model is given as follows:

$$P(Y=1\mid x_1,x_2,x_3)=\frac{{(1+e^{1.374-0.860x_1+0.355x_2-0.201x_3})}^{11.93}-1}{{(1+e^{1.374-0.860x_1+0.355x_2-0.201x_3})}^{11.93}}$$

Now, the sample is divided into two subsamples. The first one, called the training sample, corresponds to 70% of the total sample, and the second one is the prediction sample (30% of the sample). From this partition, the following results for the fitted PHLR model are obtained.

According to the results in

Table 2. Model predictive capacity.

Actual/Forecast	$\widehat{y}=1$	$\widehat{y}=0$	Total
$y=1$	57	28	85
$y=0$	0	38	38
Total	57	66	123

, the accuracy is 77.23%, the sensitivity rate is 67.05%, and the specificity rate is 100%.

On the other hand,

Table 3. Skewness and kurtosis of the PHLR model for different $\alpha$ values.

$\alpha$	0.050	0.125	0.250	0.500	0.750	1.000	1.500
Skewness	0.355	0.032	0.160	0.135	0.058	0.000	−0.081
Kurtosis	1.673	2.159	2.716	2.974	2.988	3.000	3.031
$\alpha$	2.500	5.000	10.000	20.000	30.000	50.000	100.000
Skewness	−0.179	−0.303	−0.410	−0.501	−0.546	−0.597	−0.655
Kurtosis	3.090	3.201	3.331	3.469	3.548	3.644	3.765

shows the performance of the PHLR model for different values of the $\alpha$ parameter.

Table 3 shows the skewness and kurtosis coefficients of the proportional hazard logistic model for different $\alpha$ values. The results indicate that the model can fit data with both negative and positive skewness, which is an advantage over traditional logistic models. Additionally, the PHLR model can fit data with varying degrees of kurtosis, both high and low.

Diagnostic analysis is a technique to detect possible influential observations and aberrant or extraneous data. In the case of the logistic model, this technique has certain similarities with the general diagnostic analysis of regression models. However, given that the response variable only takes the values 0 and 1, a somewhat unusual situation arises. Certain difficulties may arise if there is a large number of zeros (or ones) when one expects to find few zeros or ones, which can be a sign of a lack of fit in the model. In the case of the PHLR model, the diagnostic analysis could be carried out using the Pearson residuals,

$$\tilde{r}=\frac{y_i-{\widehat{p}}_i}{\sqrt{{\widehat{p}}_i(1-{\widehat{p}}_i)}}$$

the square of which is the ith component of the Pearson chi-square statistic, the residual deviance

$$t_{D_i}=\mathrm{sign}\left(\tilde{r}\right)\sqrt{-2\left[y_{i}log\left(\frac{y_i}{{\widehat{p}}_i}\right)+(1-y_i)log\left(\frac{1-y_i}{1-{\widehat{p}}_i}\right)\right]},$$

which is an adapted version of Cook’s distance for the case of the logistic regression model (see Christensen [25]). When $y_i=0$, $t_{D_i}=\mathrm{sign}\left(\tilde{r}\right)\sqrt{-2log\left(1-{\widehat{p}}_i\right)}$, while if $y_i=1$, $t_{D_i}=\mathrm{sign}\left(\tilde{r}\right)\sqrt{-2log\left({\widehat{p}}_i\right)}$.

For the student dropout data, the residual deviance graph for the PHLR model is presented in Figure 2. Note that in this graph, there are no observations with high values of the residuals, which indicates that the model has a good fit. Likewise, the graph of the PHL distribution for the fitted probabilities is shown. Note that there are five values falling within the +2.5/−2.5 range in Figure 2b, and six values in Figure 2c, indicating that these observations are not extremely influential. Additionally, there are no observations outside the confidence bands in the envelope graphs (Figure 3b), suggesting that the PHLR model effectively handles observations that deviate slightly from the +2/−2 range.

The $rM{T}_i$ envelope graphs generated for the logistic and proportional hazard logistic models are presented in Figure 3a and Figure 3b, respectively. It is observed that the proportional hazard logistic regression model presents a better fit than the logistic regression model.

6. Conclusions

In this work, we have proposed the PHLR, a nonlinear regression model that captures complex relationships between independent variables and the response variable, particularly in the case of dichotomous data where the relationships cannot be adequately represented by a straight line. The flexibility of the PHLR model allows for a better fit to the data compared to linear models or even the logistic model.

The information matrix of the PHLR model is non-singular, ensuring that the parameters are uniquely estimable, avoiding linear dependency among them, and allowing for the proper calculation of the variance–covariance of the estimators. This guarantees the convergence of optimization and estimation algorithms, and ensures that the maximum likelihood estimators have desirable asymptotic properties, such as asymptotic normality.

In terms of information criteria such as AIC, CAIC, and BIC, the PHLR model shows a better fit than the logistic model for the analyzed student dropout data. The logistic model is revealed as a special case of the PHLR model. Additionally, the PHLR model demonstrates a good rate of correct classifications in the studied data. An alternative for the diagnostic analysis of model errors has also been proposed, offering useful tools for its implementation in educational problems or other contexts with dichotomous responses.