Bayesian Analysis of Exponentiated-Weibull Weibull Distribution for Modelling Student Dropout Times

1. Introduction

The issue of student dropout in online education has become a critical concern for educators and institutions alike. Online education offers unparalleled flexibility and accessibility, allowing students from diverse backgrounds to pursue their academic goals. However, the flexibility that makes online learning appealing also presents unique challenges, particularly regarding student engagement and retention [1,2]. High dropout rates in online courses not only affect individual students’ academic progress but also undermine the efficacy of educational programs and the reputation of institutions offering these courses. With the rapid expansion of online learning platforms, understanding the factors that influence student retention and dropout rates has become crucial for designing effective interventions aimed at mitigating dropout rates and enhancing student success [1,2,3,4].

Traditional models of student retention, such as those used in face-to-face learning environments [5,6,7,8], often fail to capture the complex and dynamic patterns observed in online dropout data. These models typically assume linear relationships and do not account for the myriad factors that can influence a student’s decision to drop out [9]. Online learning environments introduce variables such as technological issues, varying levels of digital literacy, and the need for self-regulation and motivation, all of which can significantly impact student retention [9]. Consequently, there is a pressing need for the development of more sophisticated statistical tools that can accurately model these complexities and provide coverage for dropout time. While most of the existing studies on modelling dropout focuses of the use of machine learning classification techniques [8,10,11,12,13], only few [2] has utilized survival analysis and probability modelling to analyze students dropout time. Thus, this study introduces the Exponentiated Weibull Weibull (EWW) distribution, a five-parameter model specifically designed to better represent dropout times in online learning environments. The EWW distribution extends traditional survival analysis models by incorporating additional parameters that allow for greater flexibility in modeling the hazard function. The hazard function, which represents the instantaneous risk of dropout at any given time, can take on various shapes, including increasing, decreasing, and bathtub-shaped hazard rates. This flexibility is particularly important for capturing the uncommon patterns of student engagement and disengagement in online courses.

The motivation for modelling student dropout times by the EWW distribution comes from the requirement to include the unusual patterns of online learning engagement and disengagement. Earlier studies have applied different statistical techniques in the analysis of dropout rates, including logistic regression and Cox proportional hazards models [2,14,15,16,17]. Moreover, these models are often inadequate to properly predict the non-monotonic hazard rates of dropout that appear in educational data, where the risk of student leaving his/her course is lower at the beginning and then increases thereafter. Using the EWW distribution in this setting allows us to investigate the dropout dynamics in a more fine-grained manner. EWW distribution provides a better fit for the data by introducing parameters which is able to capture both the initial decrease, as well as long-term increase in dropout hazard rates. A closer fit is important in generating appropriate interventions that directly target dropout and retention rate for online learning courses [2].

The EWW distribution is built on the foundational principles of the Weibull distribution, a widely used model in reliability engineering and survival analysis due to its flexibility in modelling various types of data [18]. The Exponentiated Weibull (EW) distribution are known for characterization of both increasing and decreasing hazard rates, making it suitable for diverse applications [18]. The EWW distribution extends this versatility by incorporating additional parameters, allowing it to accommodate more complex hazard functions often observed in real-world data, such as the bathtub-shaped hazard function which maybe observed in dropout studies [2,19,20]. This enhancement is important for capturing the unusual patterns of student engagement and disengagement, providing a more accurate and detailed understanding of dropout dynamics in online learning environments.

Bayesian estimation has gained popularity in recent years due to its ability to incorporate prior information and provide more accurate parameter estimates, especially in small sample sizes or complex models [18,21,22,23,24,25,26]. Unlike the traditional Maximum Likelihood Estimation (MLE), which relies solely on the observed data, Bayesian methods combine prior distributions with the likelihood of the observed data to generate posterior distributions. This approach can yield more robust estimates, particularly in scenarios with high variability or limited data [27,28].

In this paper, I used both MLE and Bayesian estimation methods to estimate the parameters of the EWW distribution. The comparison of these methods is essential to determine the most effective approach for parameter estimation in this context. Simulation studies and real-life data analysis are conducted to evaluate the accuracy and precision of these estimates.

2. Exponentiated-Weibull Distribution

The exponentiated Weibull (EW) distribution was introduced by [20]. This distribution serves as an extension of the popular Weibull distribution, offering a more flexible model for analyzing lifetime data. Unlike the traditional Weibull distribution, the EW distribution can model both monotone and nonmonotone failure rates, making it particularly useful for practical applications where data often exhibit bathtub-shaped or upside-down bathtub-shaped failure rates. In many real-world scenarios, such as reliability engineering and survival analysis, lifetime data do not always follow a simple increasing or decreasing hazard function. Instead, they might show a complex pattern where the failure rate decreases initially, stabilizes, and then increases, resembling a bathtub shape. Similarly, an upside-down bathtub shape might be observed, where the failure rate initially increases, stabilizes, and then decreases. The EW distribution captures these patterns more accurately than distributions with purely monotone failure rates. The flexibility of the EW distribution is derived from its two shape parameters $\eta$ and $\lambda$ and one rate parameter $\delta$. The probability density function (pdf) of the three-parameter exponentiated Weibull distribution is given by:

$$f\left(t\right)=\eta \delta \lambda {\left(\delta t\right)}^{\eta -1}{\left(1-exp[-{\left(\delta t\right)}^{\eta }]\right)}^{\lambda -1}exp[-{\left(\delta t\right)}^{\eta }],$$

(1)

and its cumulative distribution function (CDF) is:

$$F\left(t\right)={\left(1-exp[-{\left(\delta t\right)}^{\eta }]\right)}^{\lambda },$$

(2)

where $t>0$ is within the support of the distribution, and $\eta >0$, $\delta >0$ and $\lambda >0$ are the parameters. If $\lambda =1$, the exponentiated Weibull distribution simplifies to the standard Weibull distribution, which has the pdf:

$$f\left(t\right)=\eta \delta {\left(\delta t\right)}^{\eta -1}exp[-{\left(\delta t\right)}^{\eta }],$$

(3)

and the cumulative distribution function (CDF) of the Weibull distribution is given by

$$F\left(t\right)=1-exp\left[-{\left(\delta t\right)}^{\eta }\right].$$

(4)

2.1. Moments and Median Survival Time

According to [29,30] , the rth moment of EW for positive integer $\lambda$ can be defined as:

$$E\left(X^r\right)=\lambda {\delta }^r\mathrm{\Gamma }\left(\frac{r}{\eta }+1\right)A_r\left(\eta \right)$$

where

$$A_r\left(\eta \right)=1+\sum _{j=1}^{\lambda -1}\frac{{(-1)}^j\left(\genfrac{}{}{0pt}{}{\lambda -1}{j}\right)}{{(j+1)}^{\left(\frac{r}{\eta }+1\right)}}$$

For the mean, with $r=1$,

$$E\left(X\right)=\lambda \delta \mathrm{\Gamma }\left(\frac{1}{\eta }+1\right)A_1\left(\eta \right)$$

$$E\left(X\right)=\lambda \delta \mathrm{\Gamma }\left(\frac{1}{\eta }+1\right)\left(1+\sum _{j=1}^{\lambda -1}\frac{{(-1)}^j\left(\genfrac{}{}{0pt}{}{\lambda -1}{j}\right)}{{(j+1)}^{\left(\frac{1}{\eta }+1\right)}}\right)$$

For the variance, with $r=2$,

$$E\left(X^2\right)=\lambda {\delta }^2\mathrm{\Gamma }\left(\frac{2}{\eta }+1\right)A_2\left(\eta \right)$$

$$E\left(X^2\right)=\lambda {\delta }^2\mathrm{\Gamma }\left(\frac{2}{\eta }+1\right)\left(1+\sum _{j=1}^{\lambda -1}\frac{{(-1)}^j\left(\genfrac{}{}{0pt}{}{\lambda -1}{j}\right)}{{(j+1)}^{\left(\frac{2}{\eta }+1\right)}}\right)$$

$$\mathrm{Var}\left(X\right)=\lambda {\delta }^2\mathrm{\Gamma }\left(\frac{2}{\eta }+1\right)\left(1+\sum _{j=1}^{\lambda -1}\frac{{(-1)}^j\left(\genfrac{}{}{0pt}{}{\lambda -1}{j}\right)}{{(j+1)}^{\left(\frac{2}{\eta }+1\right)}}\right)-{\left[\lambda \delta \mathrm{\Gamma }\left(\frac{1}{\eta }+1\right)\left(1+\sum _{j=1}^{\lambda -1}\frac{{(-1)}^j\left(\genfrac{}{}{0pt}{}{\lambda -1}{j}\right)}{{(j+1)}^{\left(\frac{1}{\eta }+1\right)}}\right)\right]}^2$$

Consequently, according to[31] the median survival time is derived as:

$$M\left(t\right)=\frac{1}{\delta }{\left[-log\left(1-0.5^{\frac{1}{\lambda }}\right)\right]}^{\frac{1}{\eta }}.$$

(5)

2.2. Survival, Hazard, and Cumulative Hazard Functions

The survival function, hazard function, and cumulative hazard function of the exponentiated Weibull distribution are as follows:

• Survival Function

  $$S\left(t\right)=1-{\left(1-exp\left[-{\left(\delta t\right)}^{\eta }\right]\right)}^{\lambda },$$
• Hazard Function

  $$h\left(t\right)=\frac{\eta \delta \lambda {\left(\delta t\right)}^{\eta -1}{\left(1-exp\left[-{\left(\delta t\right)}^{\eta }\right]\right)}^{\lambda -1}exp\left[-{\left(\delta t\right)}^{\eta }\right]}{1-{\left(1-exp\left[-{\left(\delta t\right)}^{\eta }\right]\right)}^{\lambda }},$$
• Cumulative Hazard Function

  $$H\left(t\right)=-log\left\{1-{\left(1-exp\left[-{\left(\delta t\right)}^{\eta }\right]\right)}^{\lambda }\right\}.$$
• Hazard Function Behavior

The behaviour of the hazard function can be characterized as follows:

• It is monotone increasing if $\eta \ge 1$ and $\eta \lambda \ge 1$.
• It is monotone decreasing if $\eta \le 1$ and $\eta \lambda \le 1$.
• It is unimodal if $\eta <1$ and $\eta \lambda >1$.
• It is bathtub-shaped if $\eta >1$ and $\eta \lambda <1$.

3. Five-Parameter Exponentiated Weibull Weibull Distribution

The cumulative distribution function (CDF) of the exponentiated Weibull-generated family [32] for a random variable $t>0$ is given by

$$F\left(t\right)={\left\{1-exp\left[-\delta {\left(\frac{G\left(t\right)}{1-G\left(t\right)}\right)}^{\eta }\right]\right\}}^{\lambda };\delta ,\eta ,\lambda >0,$$

(6)

where $\eta$ and $\lambda$ are the two shape parameters, and $\delta$ is the rate parameter. This CDF provides a broader family of continuous distributions. The corresponding probability density function (PDF) is given by

$$f\left(t\right)=\frac{\lambda \delta \eta {\left(G\left(t\right)\right)}^{\eta -1}g\left(t\right)}{{\left[1-\left(G\left(t\right)\right)\right]}^{\eta +1}}exp\left[-\delta {\left(\frac{G\left(t\right)}{1-G\left(t\right)}\right)}^{\eta }\right]{\left\{1-exp\left[-\delta {\left(\frac{G\left(t\right)}{1-G\left(t\right)}\right)}^{\eta }\right]\right\}}^{\lambda -1};\delta ,\eta ,\lambda >0.$$

(7)

In this study, we propose a novel five-parameter model that serves as a competitive extension to the Weibull distribution by utilizing the EW-G distribution. By replacing $G\left(t\right)$ in equation (6) with the cumulative distribution function of the Weibull distribution from equation (4), we obtain the following result.

$$F\left(t\right)={\left\{1-exp\left[-\delta {\left(\frac{1-exp\left[-{\left(\nu t\right)}^{\gamma }\right]}{exp\left[-{\left(\nu t\right)}^{\gamma }\right]}\right)}^{\eta }\right]\right\}}^{\lambda }$$

Simplifying the expression inside the exponential:

$$\frac{1-exp\left[-{\left(\nu t\right)}^{\gamma }\right]}{exp\left[-{\left(\nu t\right)}^{\gamma }\right]}=exp\left[{\left(\nu t\right)}^{\gamma }\right]-1$$

So the function becomes:

$$F\left(t\right)={\left\{1-exp\left[-\delta {\left(exp\left[{\left(\nu t\right)}^{\gamma }\right]-1\right)}^{\eta }\right]\right\}}^{\lambda }$$

Therefore, the cumulative distribution function (CDF) of the five-parameter Exponentiated-Weibull Weibull (EWW) distribution for a random variable $t>0$ is given by:

$$F\left(t\right)={\left\{1-exp\left[-\delta {\left(exp\left[{\left(\nu t\right)}^{\gamma }\right]-1\right)}^{\eta }\right]\right\}}^{\lambda };\delta ,\eta ,\lambda ,\nu ,\gamma >0,$$

(8)

where $\eta$, $\delta$ and $\gamma$ are the three shape parameters, and $\delta$ and $\nu$ are the rate parameters.

Correspondingly, substituting equations (3) and (4) in (7), we have the pdf of the five-parameter exponentiated-Weibull Weibull distribution given by

$$\begin{array}{cc}\hfill f\left(t\right)& =\lambda \delta \eta \gamma \nu {\left(\nu t\right)}^{\gamma -1}{\left\{1-exp\left[-{\left(\nu t\right)}^{\gamma }\right]\right\}}^{\eta -1}exp\left[\eta {\left(\nu t\right)}^{\gamma }\right]\hfill \\ \hfill & \times exp\left[-\delta {\left(exp\left[{\left(\nu t\right)}^{\gamma }\right]-1\right)}^{\eta }\right]\hfill \\ \hfill & \times {\left\{1-exp\left[-\delta {\left(exp\left[{\left(\nu t\right)}^{\gamma }\right]-1\right)}^{\eta }\right]\right\}}^{\lambda -1};\delta ,\eta ,\lambda ,\nu ,\gamma >0.\hfill \end{array}$$

(9)

The survival, hazard and cumulative hazard functions follow from the cumulative density function as follows respectively:

$$\begin{array}{cc}\hfill S\left(t\right)& =1-F\left(t\right)\hfill \\ \hfill & =1-{\left\{1-exp\left[-\delta {\left(exp\left[{\left(\nu t\right)}^{\gamma }\right]-1\right)}^{\eta }\right]\right\}}^{\lambda },\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill h\left(t\right)& =\frac{f\left(t\right)}{1-F\left(t\right)}\hfill \\ \hfill & =\frac{f\left(t\right)}{1-{\left\{1-exp\left[-\delta {\left(exp\left[{\left(\nu t\right)}^{\gamma }\right]-1\right)}^{\eta }\right]\right\}}^{\lambda }},\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill H\left(t\right)& =-log\left(1-F\left(t\right)\right)\hfill \\ \hfill & =-log\left(1-{\left\{1-exp\left[-\delta {\left(exp\left[{\left(\nu t\right)}^{\gamma }\right]-1\right)}^{\eta }\right]\right\}}^{\lambda }\right).\hfill \end{array}$$

(12)

In addition, the median survival time is given by:

$$M\left(t\right)=\frac{1}{\nu }{\left(ln\left\{1+{\left[-\frac{ln\left(1-0.5^{\frac{1}{\lambda }}\right)}{\delta }\right]}^{\frac{1}{\eta }}\right\}\right)}^{\frac{1}{\gamma }}.$$

(13)

Note that the derivations in equations (8) to (13) are different from the one presented by [32] in terms of parameterization of the Weibull base distribution. While the approach of [32] is applicable in medicine, our own parameterization is applicable for modelling failure times [33,34] which is relevant for the specific situation addressed in this paper.

Figure 1 shows that when $\lambda ,\eta ,\gamma ,\delta ,\nu$ are all equal to 1, we have an increasing hazard and when $\lambda ,\eta ,\gamma$ are less than 1, we have a decreasing hazard. Similarly, when $\lambda ,\gamma$ are less than 1 but $\gamma$ at least 1, we have a bathtub hazard and $\eta$ is less than 1 but $\lambda$ and $\gamma$ are greater than 1, we have a unimodal hazard.

4. Modelling Student Dropout Time

Suppose we have a dataset $D=\{t_i,Y_i\}$ that corresponds to dropout times $t_i$ and outcomes (withdrawn: W and completed: C) $Y_i$ for a randomly selected student i from an online class module presentation. Thus, the joint density function for simultaneously modelling dropout and completion times is given by:

$$f(t_W,t_C)={\left[f\left(W\right)\right]}^{\alpha }\times {\left[S\left(C\right)\right]}^{1-\alpha }$$

(14)

where $\alpha =\{0,1\}$ with 1 denoting withdrawn and 0 denoting complete, $f\left(W\right)$ and $S\left(C\right)$ are the EWW probability density function and survival functions for the times of withdrawn and completed students, respectively. Substituting (9) and (10) in (12) we have:

$$\begin{array}{cc}\hfill f(t_W,t_C)& =\left(\lambda \delta \eta \gamma \nu {\left(\nu t\right)}^{\gamma -1}{\left\{1-exp\left[-{\left(\nu t\right)}^{\gamma }\right]\right\}}^{\eta -1}exp\left[\eta {\left(\nu t\right)}^{\gamma }\right]\right.\hfill \\ \hfill & \times exp\left[-\delta {\left(exp\left[{\left(\nu t\right)}^{\gamma }\right]-1\right)}^{\eta }\right]\hfill \\ \hfill & {\left.\times {\left\{1-exp\left[-\delta {\left(exp\left[{\left(\nu t\right)}^{\gamma }\right]-1\right)}^{\eta }\right]\right\}}^{\lambda -1}\right)}^{\alpha }\hfill \\ \hfill & \times {\left(1-{\left\{1-exp\left[-\delta {\left(exp\left[{\left(\nu t\right)}^{\gamma }\right]-1\right)}^{\eta }\right]\right\}}^{\lambda }\right)}^{1-\alpha }.\hfill \end{array}$$

(15)

5. Parameter Estimation

5.1. Maximum Likelihood Estimation

The maximum likelihood estimates of the unknown parameter set $\mathrm{\Omega }=\{\lambda ,\eta ,\gamma ,\delta ,\nu \}$ can be obtained from the random sample of students withdraw and completion times $t_i,i\in \{1,\cdots ,n\}$. The likelihood and log-likelihood functions for the joint density in (15) are given as:

$$\begin{array}{cc}\hfill L(t_1,\cdots ,t_n;\mathrm{\Omega })& =\prod _{i=1}^{n}f(t_W,t_C;\mathrm{\Omega })\hfill \\ \hfill & =\prod _{i=1}^n\left(\lambda \delta \eta \gamma \nu {\left(\nu t\right)}^{\gamma -1}{\left\{1-exp\left[-{\left(\nu t\right)}^{\gamma }\right]\right\}}^{\eta -1}exp\left[\eta {\left(\nu t\right)}^{\gamma }\right]\right.\hfill \\ \hfill & \times exp\left[-\delta {\left(exp\left[{\left(\nu t\right)}^{\gamma }\right]-1\right)}^{\eta }\right]\hfill \\ \hfill & {\left.\times {\left\{1-exp\left[-\delta {\left(exp\left[{\left(\nu t\right)}^{\gamma }\right]-1\right)}^{\eta }\right]\right\}}^{\lambda -1}\right)}^{\alpha }\hfill \\ \hfill & \times {\left(1-{\left\{1-exp\left[-\delta {\left(exp\left[{\left(\nu t\right)}^{\gamma }\right]-1\right)}^{\eta }\right]\right\}}^{\lambda }\right)}^{1-\alpha }\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill log\left(L(t_1,\cdots ,t_n;\mathrm{\Omega })\right)& =\sum _{i=1}^{n}log\left(f(t_W,t_C;\mathrm{\Omega })\right)\hfill \\ \hfill & =\sum _{i=1}^{n}log\left[(\lambda \delta \eta \gamma \nu {\left(\nu t\right)}^{\gamma -1}{\left\{1-exp\left[-{\left(\nu t\right)}^{\gamma }\right]\right\}}^{\eta -1}exp\left[\eta {\left(\nu t\right)}^{\gamma }\right]\right.\hfill \\ \hfill & \times exp\left[-\delta {\left(exp\left[{\left(\nu t\right)}^{\gamma }\right]-1\right)}^{\eta }\right]\hfill \\ \hfill & \times {\left\{1-exp{\left[-\delta \left(exp\left[{\left(\nu t\right)}^{\gamma }\right]{-1)}^{\eta }\right]\right\}}^{\lambda -1}\right)}^{\alpha }\hfill \\ \hfill & \left.\times {\left(1-{\left\{1-exp\left[-\delta {\left(exp\left[{\left(\nu t\right)}^{\gamma }\right]-1\right)}^{\eta }\right]\right\}}^{\lambda }\right)}^{1-\alpha }\right]\hfill \end{array}$$

(17)

The MLE estimates of the parameters correspond to the solutions of the partial derivatives of the log-likelihood function concerning $\lambda$, $\eta$, $\gamma$, $\delta$, and $\nu$. Given the $log\left(L(t_1,\cdots ,t_n;\mathrm{\Omega })\right)$, we can define $f(t;\mathrm{\Omega })$ as:

$$\begin{array}{cc}\hfill f(t;\mathrm{\Omega })& =\left[(\lambda \delta \eta \gamma \nu {\left(\nu t\right)}^{\gamma -1}{\left\{1-exp\left[-{\left(\nu t\right)}^{\gamma }\right]\right\}}^{\eta -1}exp\left[\eta {\left(\nu t\right)}^{\gamma }\right]\right.\hfill \\ \hfill & \times exp\left[-\delta {\left(exp\left[{\left(\nu t\right)}^{\gamma }\right]-1\right)}^{\eta }\right]\hfill \\ \hfill & \times {\left\{1-exp{\left[-\delta \left(exp\left[{\left(\nu t\right)}^{\gamma }\right]{-1)}^{\eta }\right]\right\}}^{\lambda -1}\right)}^{\alpha }\hfill \\ \hfill & \left.\times {\left(1-{\left\{1-exp\left[-\delta {\left(exp\left[{\left(\nu t\right)}^{\gamma }\right]-1\right)}^{\eta }\right]\right\}}^{\lambda }\right)}^{1-\alpha }\right].\hfill \end{array}$$

Suppose we define $k\left(t\right)$ and $m\left(t\right)$ for simplification:

$$\begin{array}{cc}\hfill k\left(t\right)& =\lambda \delta \eta \gamma \nu {\left(\nu t\right)}^{\gamma -1}{\left\{1-exp\left[-{\left(\nu t\right)}^{\gamma }\right]\right\}}^{\eta -1}exp\left[\eta {\left(\nu t\right)}^{\gamma }\right]\hfill \\ \hfill & \times exp\left[-\delta {\left(exp\left[{\left(\nu t\right)}^{\gamma }\right]-1\right)}^{\eta }\right]\hfill \\ \hfill & \times {\left\{1-exp\left[-\delta {\left(exp\left[{\left(\nu t\right)}^{\gamma }\right]-1\right)}^{\eta }\right]\right\}}^{\lambda -1}\hfill \end{array}$$

$$m\left(t\right)=1-{\left\{1-exp\left[-\delta {\left(exp\left[{\left(\nu t\right)}^{\gamma }\right]-1\right)}^{\eta }\right]\right\}}^{\lambda }$$

Then, the log-likelihood function becomes:

$$log\left(L(t_1,\cdots ,t_n;\mathrm{\Omega })\right)=\sum _{i=1}^{n}log\left(k{\left(t_i\right)}^{\alpha }·m{\left(t_i\right)}^{1-\alpha }\right)$$

$$=\sum _{i=1}^n\left(\alpha logk\left(t_i\right)+(1-\alpha )logm\left(t_i\right)\right)$$

The corresponding partial derivatives with respect to the parameters are as follows:

1. With respect to $\lambda$:

$$\frac{\partial log\left(L\right)}{\partial \lambda }=\sum _{i=1}^n\left(\alpha \frac{1}{k\left(t_i\right)}\frac{\partial k\left(t_i\right)}{\partial \lambda }+(1-\alpha )\frac{1}{m\left(t_i\right)}\frac{\partial m\left(t_i\right)}{\partial \lambda }\right)$$

2. With respect to $\eta$:

$$\frac{\partial log\left(L\right)}{\partial \eta }=\sum _{i=1}^n\left(\alpha \frac{1}{k\left(t_i\right)}\frac{\partial k\left(t_i\right)}{\partial \eta }+(1-\alpha )\frac{1}{m\left(t_i\right)}\frac{\partial m\left(t_i\right)}{\partial \eta }\right)$$

3. With respect to $\gamma$:

$$\frac{\partial log\left(L\right)}{\partial \gamma }=\sum _{i=1}^n\left(\alpha \frac{1}{k\left(t_i\right)}\frac{\partial k\left(t_i\right)}{\partial \gamma }+(1-\alpha )\frac{1}{m\left(t_i\right)}\frac{\partial m\left(t_i\right)}{\partial \gamma }\right)$$

4. With respect to $\delta$:

$$\frac{\partial log\left(L\right)}{\partial \delta }=\sum _{i=1}^n\left(\alpha \frac{1}{k\left(t_i\right)}\frac{\partial k\left(t_i\right)}{\partial \delta }+(1-\alpha )\frac{1}{m\left(t_i\right)}\frac{\partial m\left(t_i\right)}{\partial \delta }\right)$$

5. With respect to $\nu$:

$$\frac{\partial log\left(L\right)}{\partial \nu }=\sum _{i=1}^n\left(\alpha \frac{1}{k\left(t_i\right)}\frac{\partial k\left(t_i\right)}{\partial \nu }+(1-\alpha )\frac{1}{m\left(t_i\right)}\frac{\partial m\left(t_i\right)}{\partial \nu }\right)$$

Combining these partial derivatives, we get the final expressions for the partial derivatives of the log-likelihood function with respect to $\lambda$, $\eta$, $\gamma$, $\delta$, and $\nu$. Solving these partial derivatives manually is very hard and hence in practice the Newton-Raphson procedure [34] is often used.

5.2. Bayesian Estimation

In this section, we discuss the Bayesian estimation of the parameters of the five-parameter exponentiated Weibull Weibull (EWW) distribution. We use the uninformative Uniform $(a,b)$ prior since there is no prior information or elicitation may be difficult [18]. Let the parameters of the EWW distribution be $\lambda$, $\eta$, $\gamma$, $\delta$, and $\nu$. We suggest five independent Uniform $(a,b)$ distributions for these parameters. The joint density function for the prior of these five parameters can be defined as:

$$\begin{array}{cc}\hfill f(\lambda ,\eta ,\gamma ,\delta ,\nu \mid a_1,b_1;a_2,b_2;a_3,b_3;a_4,b_4;a_5,b_5)& =\frac{1}{b_1-a_1}\times \frac{1}{b_2-a_2}\times \frac{1}{b_3-a_3}\hfill \\ \hfill & \times \frac{1}{b_4-a_4}\times \frac{1}{b_5-a_5}\hfill \end{array}$$

(18)

where $a_1,b_1;a_2,b_2;a_3,b_3;a_4,b_4;a_5,b_5$ are the prior hyperparameters for the parameters $\lambda ,\eta ,\gamma ,\delta ,\nu$ respectively. The posterior distribution of the parameters $\lambda ,\eta ,\gamma ,\delta ,\nu$ for the EWW model can be defined as the product of the likelihood $L(t_1,\cdots ,t_n;\mathrm{\Omega })$ and the prior density which is:

$$f(\lambda ,\eta ,\gamma ,\delta ,\nu \mid t,a_1,b_1;a_2,b_2;a_3,b_3;a_4,b_4;a_5,b_5)=L(t_1,\cdots ,t_n;\mid \lambda ,\eta ,\gamma ,\delta ,\nu )\times f(\lambda ,\eta ,\gamma ,\delta ,\nu )$$

(19)

where the likelihood function $L(t_1,\cdots ,t_n;\mathrm{\Omega })$ is (16).

The posterior distribution in (19) does not have a closed form as it is an approximate distribution since the marginal distribution that ensures it scales to one has been dropped. One of the ways of sampling from this distribution is by using the Metropolis-Hastings algorithm [24].

5.2.1. Metropolis-Hastings Algorithm (MHA) for EWW

The Metropolis-Hastings algorithm [24] is a hybrid random walk that employs the accept/reject algorithm to arrive at the desired distribution. The desired Bayesian inference can be summarized as follows:

Algorithm 1: Metropolis-Hastings Algorithm

6. Simulation study and results

To illustrate the EWW dropout model using both the maximum likelihood estimator and Bayesian estimator, we simulated random samples that corresponds to the four hazard types presented in Figure 1. Specifically, the following simulation parameters were used for the various hazard types.

• Monotone increasing: $\lambda =1,\eta =1,\gamma =1,\delta =1,\nu =1$.
• Monotone decreasing: $\lambda =.7,\eta =.7,\gamma =.7,\delta =1,\nu =1$.
• Unimodal: $\lambda =4,\eta =.3,\gamma =2,\delta =3,\nu =3$.
• Bathtub-shaped: $\lambda =.7,\eta =1,\gamma =.6,\delta =1,\nu =1$.

The dropout ($t_i$) time was generated using the following equation:

$$t_i=\frac{1}{\nu }{\left(ln\left\{1+{\left[-\frac{ln\left(1-U^{\frac{1}{\lambda }}\right)}{\delta }\right]}^{\frac{1}{\eta }}\right\}\right)}^{\frac{1}{\gamma }},$$

(20)

where U is the standardized uniform random variable with parameters $0,1$. To achieve censoring, a random censoring scheme was adopted. In this scheme, a random variable Z that follows a Bernoulli distribution with parameter $\pi$, denoted as Bernoulli$\left(\pi \right)$, is defined. If $Z_i=1$, then $t_i$ corresponds to the dropout time of student i. Otherwise, $t_i$ corresponds to the completion time of student i. To study the effect of censoring on the estimators, we set $\pi =0,0.1,0.2,0.3,0.4,0.5$ which corresponds to $10\%-50\%$ censoring percentages. In addition, the estimators were evaluated using bias and root mean square error. Bias and root mean square error for each parameter set are given below:

$$Bias\left(\theta \right)=\widehat{\theta }-\theta ,$$

(21)

$$Bias\left(\mathrm{\Omega }\right)=\frac{{\sum }_{i=1}^5(\widehat{\mathrm{\Omega }}-\mathrm{\Omega })}{5},$$

(22)

$$RMSE\left(\theta \right)=\sqrt{{\left(\widehat{\theta }-\theta \right)}^2},$$

(23)

$$RMSE\left(\mathrm{\Omega }\right)=\frac{\sqrt{{\sum }_{i=1}^5{\left(\widehat{\mathrm{\Omega }}-\mathrm{\Omega }\right)}^2}}{5}.$$

(24)

The sample size n was fixed at 200 throughout the experiment. To achieve stability, each simulation run was replicated 1000 times [35]. For the Bayesian estimation, $20,000$ posterior samples were generated with a thinning value set to 10 (i.e., every 10th estimate was saved). The initial burn-in period, during which samples were discarded, was set at 2000.

In Table 1, we present the parameter estimates, standard deviations (SD) for Bayesian estimates, standard errors (SE) for Maximum Likelihood Estimates (MLE), biases, and root mean squared errors (RMSE) for both Bayesian and MLE methods. The table summarizes these statistics for five different parameters: $\lambda$, $\eta$, $\gamma$, $\delta$, and $\nu$, evaluated at different censoring levels ($\pi$) of 10%, 20%, 30%, 40%, and 50% for scenario 1: increasing hazard type. The results show that Bayesian estimates consistently outperform Maximum Likelihood Estimates (MLE) across all parameters. For $\lambda$, Bayesian estimates are close to 1 the true value with low standard deviations and root mean squared errors (RMSE), whereas MLEs exhibit higher variability and larger biases, particularly at the 10% quantile. For $\eta$, Bayesian estimates also display consistency and low biases, while MLEs vary significantly with large biases and higher RMSE values.

Similarly Table 1 shows that the Bayesian estimates for $\gamma$ are accurate with small positive biases and low RMSE values, whereas MLEs show greater variability and higher biases, notably at the 20% quantile. For $\delta$, Bayesian estimates are close to the true value with low biases and RMSE, while MLEs are more variable with larger biases and higher RMSE values. Lastly, Bayesian estimates for $\nu$ are accurate with small biases and low RMSE, in contrast to MLEs which exhibit higher variability, larger biases, and higher RMSE values. Overall, Bayesian estimation demonstrates higher accuracy and precision compared to MLE for all parameters, indicating its robustness in this context.

Table 2 presents the parameter estimates, standard deviations (SD) for Bayesian estimates, standard errors (SE) for Maximum Likelihood Estimates (MLE), biases, and root mean squared errors (RMSE) for both Bayesian and MLE methods evaluated at different censoring levels ($\pi$) of 10%, 20%, 30%, 40%, and 50% for a decreasing hazard type. The results show that Bayesian estimates consistently exhibit lower biases and RMSE values compared to MLEs across all parameters and censoring levels. Bayesian estimates for $\lambda$, $\eta$, $\gamma$, $\delta$, and $\nu$ are more accurate and precise, with smaller SDs and biases. In contrast, MLEs show higher variability, larger biases, and greater RMSE values, indicating less reliable estimates. For instance, the bias for $\lambda$ in Bayesian estimates remains small and close to zero, while MLEs show a significant bias, especially at the 10% censoring level. Similarly, Bayesian estimates for $\eta$ and $\nu$ demonstrate consistent accuracy, whereas MLEs vary widely with substantial biases and higher RMSE values.

Table 3 presents parameter estimates, standard deviations (SD) for Bayesian estimates, standard errors (SE) for Maximum Likelihood Estimates (MLE), biases, and root mean squared errors (RMSE) for both Bayesian and MLE methods, evaluated at different censoring levels ($\pi$) of 10%, 20%, 30%, 40%, and 50% for a unimodal hazard type. For parameter $\lambda$, Bayesian estimates are consistently close to 4 with low SD and RMSE, whereas MLEs show significant variability and larger biases, especially at the 10% censoring level. Bayesian estimates for $\eta$ remain close to the true value of 0.3, displaying small biases and low RMSEs, while MLEs vary considerably with large biases and higher RMSEs. For $\gamma$, Bayesian estimates are accurate with low biases and RMSEs, but MLEs exhibit higher variability and significant biases. In the case of $\delta$, Bayesian estimates are near the true value of 3 with low biases and RMSEs, while MLEs show greater variability and higher biases. Lastly, Bayesian estimates for $\nu$ demonstrate accuracy with small biases and low RMSEs, in contrast to MLEs which display high variability, larger biases, and elevated RMSEs.

The results presented in Table 4 indicate that Bayesian estimates consistently outperform Maximum Likelihood Estimates (MLE) across various parameters and censoring levels for the bathtub hazard type. Bayesian estimates for all parameters ($\lambda$, $\eta$, $\gamma$, $\delta$, $\nu$) show low standard deviations (SD) and root mean squared errors (RMSE), with biases generally close to zero. In contrast, MLEs exhibit higher variability, larger biases, and higher RMSE values across the board. For example, Bayesian estimates of $\lambda$ at the 10% censoring level have a low bias of 0.024 and an RMSE of 0.029, whereas MLEs show a high bias of 0.152 and an RMSE of 0.729. Similarly, Bayesian estimates of $\delta$ are close to the true value with low biases and RMSEs across all censoring levels, while MLEs are more variable with larger biases and higher RMSE values.

Overall across the four hazard types and censoring levels, the Bayesian estimates demonstrate higher accuracy and precision compared to the MLEs for all parameters, as indicated by lower biases and RMSE values. The standard deviations for Bayesian estimates are also consistently smaller than the standard errors for MLEs, further highlighting the robustness of the Bayesian approach in this context. This simulation underscores the advantages of Bayesian estimation in providing more reliable parameter estimates for the Exponentiated Weibull Weibull (EWW) distribution.

7. Application to Modelling Time to Students’ Dropout from an Online Class Presentation

The dataset used in this study was extracted from the student behaviour and performance dataset for 23,344 students in the United Kingdom, accessible at [36]. This dataset includes various data points such as assessment results and extensive logs of student interactions with the Virtual Learning Environment (VLE). The final result of each student is categorized as Distinction, Pass, Withdrawn, or Fail, along with the submission date measured in days since the start of the module presentation. For the analysis in this study, a strategic sample of 500 students was selected to reflect the actual distribution of the final results. The goodness of fit for the proposed EWW distribution, alongside the Exponential Weibull (EW), Weibull Weibull (WW), and Weibull Exponential (WE) distributions [32], was assessed using the log-likelihood, Akaike Information Criterion (AIC), and Bayesian Information Criterion (BIC). These three competing distributions were selected because they are special cases of the EWW distribution. The WW can be derived from the EWW if we set $\lambda =1$. The EWW reduces to the WE if $\lambda =1$ and $\gamma =1$. Similarly, the EWW reduces to the EW if $\lambda =1$ and $\eta =1$.

Figure 2 shows the plots of $cdf$, $pdf$, for the five-parameter exponentiated-Weibull Weibull (EWW) alongside the Exponential Weibull (EW), Weibull Weibull (WW), and Weibull Exponential (WE) distributions for the student dropout dataset based on the Bayesian parameter estimates. The plot shows that the most appropriate distribution for the student dropout dataset is EWW.

The results presented in Table 5 highlight the superior performance of Bayesian estimates over Maximum Likelihood Estimates (MLE) for various parameters for the student dropout dataset. The Bayesian estimates for $\lambda$, $\eta$, $\gamma$, $\delta$, and $\nu$ exhibit lower standard deviations (SD), indicating higher precision. For instance, the Bayesian estimate of $\lambda$ is 0.233 with an SD of 0.008, compared to the MLE estimate of 0.536 with a standard error (SE) of 0.137. Similarly, Bayesian estimates for $\eta$ and $\gamma$ are close to their true values with minimal SDs (0.002), while MLE estimates show greater variability and higher SEs (0.149 and 0.154, respectively). For the parameter $\delta$, the Bayesian method yields an estimate of 0.635 with an SD of 0.005, significantly more precise than the MLE estimate of 0.759 with an SE of 0.129. Lastly, the Bayesian estimate for $\nu$ is 0.358 with an SD of 0.003, whereas the MLE estimate is 1.190 with a considerably higher SE of 0.196. Overall, these findings underscore the robustness and precision of Bayesian estimation in parameter estimation compared to the MLE method for the student dropout dataset.

The goodness-of-fit results presented in Table 6 for the five-parameter Exponentiated Weibull Weibull (EWW) model, alongside the Exponential Weibull (EW), Weibull Weibull (WW), and Weibull Exponential (WE) models for the student dropout dataset, indicate that the EWW model provides the best fit. This is evidenced by the EWW model having the highest LogLikelihood (-750.75) and the lowest values for the Akaike Information Criterion (AIC = 1511.5) and the Bayesian Information Criterion (BIC = 1532.57). In contrast, the WW, WE, and EW models exhibit significantly lower LogLikelihood values (-1799.22, -1859.54, and -1029.02 respectively) and higher AIC (3606.44, 3725.09, and 2064.04 respectively) and BIC (3623.3, 3737.73, and 2076.68 respectively) values. These findings underscore the superior performance of the EWW model in capturing the underlying characteristics of the student dropout data, making it the most suitable choice among the compared models.

Table 7 shows the estimate of the dropout times for the student from the online class module using the four distributions based on the Bayesian parameter estimate. The result showed that median dropout time using the EWW distribution from the online class presentation is 0.8 months (24 days) from the day the class commenced. On the other hand, if any of the WW, WE, EW distributions is assumed, the median dropout time from the online class presentation is approximately 2 months (60 days) from the day the class commenced. The median estimate using any of WW, WE, EW distribution is significantly larger than EWW and also unrealistic given that most students dropout from online class modules in the first few days after the commencement of the online presentation [11,15].

8. Conclusion

In this paper, the five-parameter Exponentiated Weibull Weibull (EWW) distribution with Weibull base failure rate parameterization was presented. Specifically, the EWW distribution was developed to model the dropout times of students in an online class module. The parameters of the model were estimated using both Maximum Likelihood and Bayesian estimation procedures to recommend the best estimation method. The simulation results revealed that, at various hazard types and censoring levels, the Bayesian estimates are more accurate and precise than the Maximum Likelihood estimates. Furthermore, the real-life analysis corroborated the findings from the simulation study, with the standard deviation of the Bayesian estimates found to be significantly lower than that of the Maximum Likelihood estimates, indicating that the Bayesian estimates are more efficient for the EWW distribution. Additionally, the comparative goodness-of-fit of the EWW distribution alongside Weibull Weibull (WW), Weibull Exponential (WE), and Exponential Weibull (EW) distributions was examined. The goodness-of-fit results showed that the EWW distribution fits the student dropout times dataset better than the three competing nested distributions. Finally, the estimate of the median dropout time of the students using the EWW distribution revealed that the median dropout time from the online presentation is 24 days. This indicates that the instructor should pay more attention to student attendance around this time.