J. Nig. Soc. Phys. Sci. 5 (2023) 871 Journal of the Nigerian Society of Physical Sciences On Bivariate Nadarajah-Haghighi Distribution derived from Farlie-Gumbel-Morgenstern Copula in the Presence of Covariates Yakubu Aliyua,∗, Umar Usmanb aDepartment of Statistics, Ahmadu Bello University Zaria-Nigeria bDepartment of Mathematics, Usmanu Danfodiyo University Sokoto-Nigeria Abstract An important alternative distribution to the Weibull, generalized exponential and gamma distributions that is used in survival analysis is the Nadarajah-Haghighi exponential distribution. Similar to the Weibull, generalized exponential and gamma distributions, the Nadarajah-Haghighi exponential distribution is an extension of the well known exponential distribution. In this paper, a copula function commonly used to model very weak linear dependence was used to introduced a bivariate Nadarajah-Haghighi distribution. The joint survival function, joint probability density function and joint cumulative distribution were given in closed form. Bayesian method of estimation was used to estimate the model parameters considering the presence of right censoring and covariates. Posterior summaries of interest were obtained via standard Markov Monte Carlo (MCMC) technique. Two real data sets were used to illustrate the importance and flexibility of the bivariate model in comparison with some competing models. It was observed that, the bivariate Nadarajah-Haghighi distribution provides a better flt than bivariate exponential, bivariate Weibull, bivariate generalized exponential and bivariate modified Weibull distributions. DOI:10.46481/jnsps.2023.871 Keywords: Exponential Distribution, Nadarajah-Haghighi Distribution, Bivariate Models, Copula Function Article History : Received: 17 June 2022 Received in revised form: 10 August 2022 Accepted for publication: 23 August 2023 Published: 28 April 2023 c© 2023 The Author(s). Published by the Nigerian Society of Physical Sciences under the terms of the Creative Commons Attribution 4.0 International license (https://creativecommons.org/licenses/by/4.0). Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI. Communicated by: O. Adeyeye 1. Introduction Exponential distribution is a well known distribution due to the constant hazard rate and memory less property it exhibits. However, in survival/ reliability studies, choosing the exponential distribution may be inappropriate since its hazard rate does not show monotone and/ or non-monotone failure rate behaviours [1]. To solve this problem, many ∗Corresponding author tel. no: +2347062898441 Email address: aliyuyakubu40@mail.com (Yakubu Aliyu ) generalizations of the exponential distribution have been developed by researchers so as to add some flexibility to the exponential distribution. These include the generalized exponential distribution developed by [2], Beta-exponential by [3], Nadarajah-Haghighi by [4]. These generalizations are among the generalizations that have received the most attention in the literature as compared with other extensions of the exponential distribution. Other generalizations of the exponential distribution include the Weibull-Burr III by [5], the generalized modified Weibull by [6], log-beta Weibull by [7], two parameter Burr X by [8] and Weibull Kumaraswamy 1 Aliyu & Usman / J. Nig. Soc. Phys. Sci. 5 (2023) 871 2 distribution by [9] to mention but a few. An extension of the exponential distribution which serves as an alternative to the Weibull, generalized exponential and gamma distributions called the Nadarajah-Haghighi exponen- tial distribution was introduced by [4]. The probability den- sity function (pdf ) and cumulative distribution function (cdf ) of the univariate Nadarajah-Haghighi exponential (NH) distri- bution with parameters θ and φ are respectively given by: f (t/θ,φ) = θφ (1 + θt)φ−1 ex p ( 1 − (1 + θt)φ ) (1) and F (t/θ,φ) = 1 − ex p ( 1 − (1 + θt)φ ) (2) where θ,φ > 0 are scale and shape parameters respectively. The pdf in (1) reduces to the exponential distribution when φ = 1. The shape of the NH density can be decreasing and unimodal, while that of the hazard rate function can be monotonically in- creasing, monotonically decreasing or constant. The survival and hazard rate functions of the NH distribution are respectively given by: S (t/θ,φ) = ex p ( 1 − (1 + θt)φ ) (3) and h (t/θ,φ) = θφ (1 + θt)φ−1 (4) The distribution has been extended by different researchers in different directions. For instance, [10] extended the distribution to the unit NH distribution which could be used to model effectively data related to rates and proportions with excess of ones. [11] extended it to the Poisson NH distribution that can be used to model reliability systems and [11] explored different classical methods of estimation to estimate the parameters of the Poisson NH distribution. [12] extended it to a three param- eter discrete NH distribution that could serves as an alternative to the Poisson, negative binomial and zero inflated Poisson distribution. It was also extended to the Nadarajah-Haghighi mixture cure rate model by [13], a model that could be used to model effectively information from a population of a mixture of two type of individuals: susceptible and cured individuals. However, in a situation where we are studying two life- times T1 and T2 associated to each unit (individual), these extensions could not be used. To solve this problem, some existing bivariate distributions in the literature such as bi- variate Weibull by [14, 15, 16, 17, 18], bivariate modified weibull by [19], bivariate exponentiated discrete Weibull distribution by [20], bivariate generalized exponential by [21, 22, 23, 24, 25, 26], bivariate generalized Rayleigh by [27], bivariate inverse weibull distribution by [28, 29, 30, 31, 32, 33], bivariate frechet by [34], bivariate inverted Topp-Leone distri- bution by [35], bivariate discrete NH by [36] and bivariate NH by [37] distributions could be used. The bivariate discrete NH distribution by [36] could only be used in modelling bivariate discrete lifetimes while the bivariate NH distribution by [37] could be used to modeled bivariate lifetimes with positive and negative dependence structure. In this paper, we generalized the exponential distribution to the bivariate exponential distribution considering the Nadarajah-Haghighi exponential distribution using the Farlie - Gumbel - Morgestern copula function. Hence, the paper aimed to introduced a bivariate Nadarajah-Haghighi exponential distribution that could be used for describing bivariate data that have weak correlation between variables in lifetime data. The distribution will extends the univariate exponential and Nadarajah-Haghighi exponential distributions, and serves as a good alternative to several bivariate distributions such as: bivariate exponential, bivariate generalized exponential and bivariate weibull distributions for modelling real-valued data. An important motivation of the article is to develop a guideline for estimating the parameters of the introduced distribution in the presence of censoring and covariates, which may be of deep interest to statisticians and practitioners. Some mathematical and statistical properties of the distribution was discussed. Bayesian method of estimation was employed in estimating the parameters of the distribution and finally, we evaluate the performance of the estimators by analyzing some real life data sets. The rest of the paper is organized as follows: in section 2, we derive the survival function, the probability density function and cumulative distribution function of the bivariate Nadarajah-Haghighi distribution, also in this section, param- eters of the distribution were estimated using the Bayesian method of estimation procedure. Application of the introduced methodology is given in section 3 and we finally conclude in section 4. 2. Bivariate Nadarajah-Haghighi Distribution In this section, the bivariate Nadarajah-Haghighi distribu- tion was developed and study some of its statistical properties. 2.1. The Model Copula functions are used in connecting the joint distribu- tion function of two or more univariate distributions. The cop- ula function is said to be bivariate when it connects the joint distribution function of only two univariate distributions. Let S (tk) be the univariate survival function for the random vari- able Tk, k = 1, 2, the joint survival function S (t1, t2) is defined as: S (t1, t2) = Cλ (S (t1) S (t2)) (5) for t1, t2 > 0 where λ is a measure of dependence between the random variables T1 and T2, and C is a copula function. The Farlie-Gumbel-Morgenstern copula (FGM) was first proposed by [38] and later by [39] and [40]. The joint survival function considering the FGM copula for random variables T1 and T2 is given by: S (t1, t2) = S (t1) S (t2) [1 + λ (1 − S (t1)) (1 − S (t2))] (6) where the dependence parameter lies between ±1 inclusive. The joint survival function in (6) reduces to the survival 2 Aliyu & Usman / J. Nig. Soc. Phys. Sci. 5 (2023) 871 3 function of the product copula function that is S (t1, t2) = S (t1) S (t2) when the dependence parameter takes the value of zero. In this case, T1 and T2 are said to be independent. The dependence parameter λ is related to the kendall and spearman rank correlation coefficients by: τ (λ) = 2λ 9 and ρ (λ) = λ 3 respectively. Hence, the FGM copula is only appropriate in modelling weak dependences. The FGM copula is useful especially when dependence between the two marginal is modest in magnitude [41]. Assume T1 and T2 be two lifetimes associated to the same individual/ device with a dependence structure given by FGM copula function. The density functions of the marginal distributions for the lifetimes T1 and T2 are given by: f1 (t1) = θ1φ1 (1 + θ1t1) φ1−1 ex p ( 1 − (1 + θ1t1) φ1 ) (7) and f2 (t2) = θ2φ2 (1 + θ2t2) φ2−1 ex p ( 1 − (1 + θ2t2) φ2 ) (8) respectively, while the survival functions are given by: S 1 (t1) = ex p ( 1 − (1 + θ1t1) φ1 ) (9) and S 2 (t2) = ex p ( 1 − (1 + θ2t2) φ2 ) (10) respectively. The joint survival function based on the FGM cop- ula given by (6) using the marginal survival functions in (9) and (10) is given by: S (t1, t2) = ex p ( 2 − (1 + θ1t1) φ1 − (1 + θ2t2) φ2 ) [ 1 + λ ( 1 − ex p ( 1 − (1 + θ1t1) φ1 )) ( 1 − ex p ( 1 − (1 + θ2t2) φ2 ))] (11) To obtain the joint density function for the random variables T1 and T2, we obtain the second derivative of S (t1, t2) with respect to t1 and t2. That is, f (t1, t2) = ∂2 S (t1,t2 ) ∂t1∂t2 which yields: f (t1, t2) = θ1θ2φ1φ2 pq1q2 [ 1 + λ (1 + 4 p − 2p1 − 2p2) ] (12) where p = ex p ( 2 − (1 + θ1t1) φ1 − (1 + θ2t2) φ2 ) , p1 = ex p ( 1 − (1 + θ1t1) φ1 ) , p2 = ex p ( 1 − (1 + θ2t2) φ2 ) , q1 = (1 + θ1t1) φ1−1 and q2 = (1 + θ2t2) φ2−1 The joint cumulative distribution function is easily obtain as: F (t1, t2) = p (1 + λ(1 − p1)(1 − p2)) − p1 − p2 + 1 (13) The marginal distribution functions are given by F (t1) = 1 − S (t1) and F (t2) = 1 − S (t2). The first partial derivatives with respect to t1i and t2i are given by: ∂S (t1i, t2i) ∂t1i = −θ1φ1 p1 p [ 1 + λ (1 − p1) (1 − p2) ] + λθ1φ1 pq1 p1 (1 − p2) (14) and ∂S (t1i, t2i) ∂t2i = −θ2φ2q2 p [ 1 + λ (1 − p1) (1 − p2) ] + p (1 − p1) p2 (15) respectively. The survival function in (11) reduce to the survival function of the product bivariate Nadarajah-Haghighi distribu- tion when the dependent parameter assumed the value zero. It also reduce to the FGM bivariate exponential distribution when φ1 = φ2 = 1 and to the product bivariate exponential distribu- tion when φ1 = φ2 = 1 and λ = 0. 2.2. Estimation In this section, the problem of estimating the parameters of the Farlie-Gumbel-Mogenstern bivariate Nadarajah-Haghighi (FGMBNH) distribution based on random samples of size n was addressed using the Bayesian estimation procedure. Let (T11, T21), (T12, T22) , · · · , (T1n, T2n) be bi- variate random sample of size n from the FGMBNH distribution, let w = (θ1, θ2, φ1, φ2, λ)′ be the vector of parameters. Then, the likelihood func- tion L (w) when the lifetimes (T11, T21), (T12, T22) , · · · , (T1n, T2n) is assumed to be non-censored can be expressed as: L (w) = n∏ i=1 f (t1i, t2i) (16) substituting equation (12) and taking natural logarithm gives: `(Θ) = nlog(θ1) + nlog(θ2) + nlog(φ1) + nlog(φ2)+ 2n − n∑ i=1 [ (1 + θ1t1i) φ1 + (1 + θ2t2i) φ2 ] + (φ1 − 1) n∑ i=1 log (1 + θ1t1i) + (φ2 − 1) n∑ i=1 log (1 + θ2t2i) + n∑ i=1 [ 1 + λ [ 1 + 4ex p ( 2 − (1 + θ1t1i) φ1 − (1 + θ2t2i) φ2 ) − 2ex p ( 1 − (1 + θ1t1i) φ1 ) − 2ex p ( 1 − (1 + θ2t2i) φ2 )]] (17) On the other hand, assume the lifetimes T1 or T2 or both T1 and T2 may be right censored. Assume also that, the cen- soring is independent of the time to the event of interest in the study. Let (T11, T21) , (T12, T22) , · · · , (T1n, T2n) be a ran- dom sample from the FGMBNH distribution with parameter w where w = (θ1,θ2,φ1,φ2,λ)′ is a parameter space. Then, each ith observation i = 1, 2, · · · , n fall in one of the following groups: • U1 : both t1i and t2i are uncensored observations. • U2 : t1i is uncensored and t2i is censored observation. • U3 : t1i is censored and t2i is uncensored observation. • U4 : both t1i and t2i are censored observations. 3 Aliyu & Usman / J. Nig. Soc. Phys. Sci. 5 (2023) 871 4 Then, the likelihood based on these conditions can be expressed as: L = ∏ i∈U1 [ ∂2S (t1i, t2i) ∂t1i∂t2i ] ∏ i∈U2 [ −∂S (t1i, t2i) ∂t1i ] × ∏ i∈U3 [ ∂S (t1i, t2i) ∂t2i ] ∏ i∈U4 S (t1i, t2i) (18) Assume τ1i and τ2i be indicator variables, such that{ τki = 1 when tki is an uncensored observation τki = 0 when tki is censored observation Then, the likelihood function in equation (18) is written as: L = n∏ i=1 [ ∂2S (t1i, t2i) ∂t1i∂t2i ]τ1iτ2i [ −∂S (t1i, t2i) ∂t1i ]τ1i (1−τ2i ) × [ −∂S (t1i, t2i) ∂t2i ](1−τ1i )τ2i [S (t1i, t2i)] (1−τ1i )(1−τ2i ) (19) observe that, the likelihood function in expression (19) reduces to the likelihood function of the non-censored observation in expression (16) when τ1i = τ2i = 1. Furthermore, in the presence of m covariates x1, x2, · · · , xm affecting the parameters of the FGMBNH distribution, a link function is assumed for the parameters φ1,φ2,θ1 and θ2. That is, φ1 = ex p (φ10 + φ11 x1 + φ12 x2 + · · · + φ1m xm) (20) φ2 = ex p (φ20 + φ21 x1 + φk2 x2 + · · · + φ2m xm) (21) θ1 = ex p (θ10 + θ11 x1 + θ12 x2 + · · · + θ1m xm) (22) and θ2 = ex p (θ20 + θ21 x1 + θ22 x2 + · · · + θ2m xm) (23) Therefore, to obtained inferences from the FGMBNH distri- bution in the presence of covariate(s), values of φ1, φ2, θ1 and θ2 in equations (20), (21), (22) and (23) respectively are appropriately substituted in the model. In Bayesian method of estimation, the joint posterior dis- tribution of the model parameters is proportional to the product of the joint prior distribution of the parameters and the likelihood function for w given by equation (16) when the lifetimes T1i and T2i are uncensored and by equation (19) when the lifetimes T1i and T2i are censored. Thus, assume the prior distributions π11(θ1) ∝ θ b1−1 1 e a1θ1 θ1 > 0 π12(θ2) ∝ θ b2−1 2 e a2θ2 θ2 > 0 π21(φ1) ∝ φ d1−1 1 e c1φ1 φ1 > 0 and π22(φ2) ∝ φ d2−1 2 e c2φ2 φ2 > 0 for the parameters θ1, θ2,φ1 and φ2 respectively, while the prior distribution of the dependence parameter (λ) was assumed to be uniform distribution. That is π3(λ) ∼ Uni f ( p, q), where a1, a2, b1, b2, c1, c2, d1, d2, p and q are known hyper- parameters. The hyper-parameters a1, a2, b1 b2, c1, c2, d1 and d2 are assumed to be non-negative while the parameters p and q are assumed to be between ±1. Lets further assumed prior independence and let the joint prior distribution for the parameters θ1, θ2,φ1, φ2 and λ be Π(w) = π11(θ1)π12(θ2)π21(φ1)π22(φ2)π3(λ) (24) Moreover, assumed normal prior distribution for the covariate parameters. That is φ jk ∼ N(µ,σ2) and θ jk ∼ N(µ,σ2) for j = 1, 2 and k = 1, 2, · · ·m, where µ and σ2 are known hyper- parameters. Hence, the joint posterior density function for w is given by: `(w/x) ∝ L × Π(w) (25) where x = (t1, t2,τ1,τ2)′ is the vector of observed lifetimes, L is the likelihood function given by equations (16) and (19), while Π(w) is the product of the product of the prior probability density functions. The full conditional posterior distributions are obtained by deriving the posterior distribution of each parameter given the data and all other parameters of the model. Posterior summaries of interest were obtained by using Markov Chain Monte Carlo (MCMC) technique. In our applications, 220,000 Gibbs samples for each model parameter was generated. However, to minimized the effect of initial values, the first 20,000 simulated samples were discarded as burn-in. Furthermore, each 20th simulated sample was stored so as to avoid auto-correlation between successive samples. The Bayesian estimates for the parameters were obtained using the medians of the respective posterior distributions since some simulated distributions were quite skewed. Credible intervals were also determine for each model parameter from the 2.5th and 97.5th centiles of the posterior distribution of each model parameter. Auto-correlation and trace plots were used in assessing the convergence of simulated samples. Different formulations were assessed using the Deviance Information Criteria (DIC). The DIC a generalization of the Akaike Information Criteria for the Bayesian anal- ysis, is obtained from the samples generated from the MCMC simulation [42]. According to [43], the DIC is computed as: DIC = D(ŵ) + 2np = 2D̄ − D(ŵ), where D(ŵ) is the deviance evaluated using the mean of the model parameters which is obtained from the MCMC samples, D̄ is the posterior mean of the deviance and np is the effective number of the model parameters which is computed as np = D̄ − D(ŵ). Better model fits are indicated by lower DIC values. All Bayesian parameter estimates, their 95% credible Interval (CrI) and the DIC for each formulation were obtained using the OpenBUGS software version 3.2.3. 4 Aliyu & Usman / J. Nig. Soc. Phys. Sci. 5 (2023) 871 5 3. Applications In this section, infections in kidney patients data from [44] and Tobacco-stained-fingers data set from [45] were used in demonstrating the applicability of the introduced methodology. The FGMBNH model is compared with its special case (prod- uct bivariate Nadarajah-Haghighi (PBNH) model) and that of some competing bivariate distributions. 3.1. Kidney data The kidney data show the recurrence times to infection at point of insertion of catheter using portable dialysis equipment. Two recurrence times were recorded for each patient together with censoring indicator (Infection occurs =1 and censored=0) and the risk variable values: age, sex (male=1, female=2) and disease type. Assume T1 and T2 refers to first and second recurrence time respectively. The data was first fitted to the FGMBNH distribution in the presence of right censoring not in- cluding covariate and compared its performance with the fits of FGM bivariate exponential (FGMBE), FGM bivariate Weibull (FGMBW), FGM bivariate generalized exponential (FGM- BGE) and FGM bivariate modified weibull (FGMBMW) distributions. It is also fitted to the special case of these bivariate distributions. That is, the prod- uct bivariate Nadarajah-Haghighi (PBNH), product bi- variate exponential (PBE), product bivariate Weibull (PBW), product bivariate generalized exponential (PBGE) and product bivariate modified Weibull (PBMW) distributions. In addition, the data is then fitted to the FGMBNH model in the presence of right censoring and covariates. In analyzing this data set, the procedure discussed in section 2.2 was followed. To be specific with the prior distributions, we assumed Gamma(1, 1) for θ1,θ2,φ1 and φ2 while we assume U(−1, 1) for the dependence parameter(λ). Furthermore, N(0, 10) was assumed for the covariate parameters. Table 1 give the posterior summary statistics for the bi- variate Nadarajah-Haghighi distribution considering the FGM copula function compared to the aforementioned bivariate distributions. The table also gives the MC error for each posterior estimate. The MC error of the posterior estimates of the fitted models were all less than 120 of standard deviation of the estimates, this showed that the posterior estimates have reasonably good precision. The results also showed that, the FGMBNH distribution is more efficient than the FGMBE, FGMBW, FGMBGE and FGMBMW distributions, since it has the least information criteria value. Furthermore, the estimates of the dependent parameter (λ) are similar for all the distributions considered. Table 2 give the posterior summary statistics for the bivari- ate Nadarajah-Haghighi distribution considering the product copula function compared to the aforementioned bivariate distributions. Similar to the results for the distributions based on the FGM copula, the FGMBNH distribution is more efficient than the aforementioned distributions since it has the least information criteria value. Comparing between the estimates of the bivariate distributions considering the FGM and product Table 1. posterior summaries considering FGM copula - Kidney data model para- med- sd MC 95% DIC meter ian error CrI FGM λ 0.5556 0.3546 0.0059 (-0.3138, 0.9781) 682.5 BNH φ1 0.4526 0.1496 0.0039 (0.2727, 0.8525) φ2 0.6546 0.2932 0.0079 (0.3430, 1.5020) θ1 0.0351 0.0330 0.0007 (0.0094, 0.1226) θ2 0.0152 0.0139 0.0003 (0.0039, 0.0535) FGM λ 0.5225 0.3345 0.0067 (-0.2720,0.9676) 686.5 BE θ1 0.0077 0.0013 0.0000 (0.0053,0.0106) θ2 0.0074 0.0015 0.0000 (0.0050,0.0106) FGM β1 0.7454 0.1022 0.0024 (0.5560,0.9589) 686.4 BW β2 0.8874 0.1225 0.0028 (0.6599,1.1390) λ 0.5608 0.3536 0.0073 (-0.3120,0.9780) θ1 0.0292 0.0190 0.0005 (0.0090,0.0797) θ2 0.0134 0.0108 0.0002 (0.0035,0.0433) FGM β1 0.7479 0.1595 0.0047 (0.4950,1.1240) 687.8 BGE β2 1.0470 0.2368 0.0072 (0.6696,1.5840) λ 0.5318 0.346 0.0105 (-0.2887,0.9743) θ1 0.0062 0.0016 0.0000 (0.0036,0.0098) θ2 0.0079 0.0020 0.0000 (0.0045,0.0122) FGM α1 0.0427 0.0211 0.0016 (0.0134,0.0942) 689.3 BMW α2 0.0217 0.0189 0.0013 (0.0047,0.0770) β1 0.6246 0.1068 0.0084 (0.4233,0.8568) β2 0.7354 0.1594 0.0128 (0.4085,1.0780) λ 0.0009 0.0007 0.0000 (0.0001,0.0026) φ1 0.0011 0.0009 0.0001 (0.0001,0.0035) φ2 0.6511 0.3084 0.0263 (-0.1261,0.9863) Table 2. posterior summaries considering product copula - Kidney data model para- med- sd MC 95% DIC meter ian error CrI PB φ1 0.4545 0.1598 0.0043 (0.2727, 0.8799) 683.2 NH φ2 0.6627 0.2860 0.0072 (0.3551, 1.4520) θ1 0.0346 0.0308 0.0005 (0.0092, 0.1237) θ2 0.0149 0.0124 0.0002 (0.0041, 0.0513) PBE θ1 0.0077 0.0013 0.0000 (0.0054,0.0106) 687.4 θ2 0.0076 0.0015 0.0000 (0.0051,0.0110) PBW β1 0.7487 0.0972 0.0034 (0.5675,0.9525) 686.9 β2 0.8998 0.1298 0.0047 (0.6717,1.1770) θ1 0.0286 0.0172 0.0005 (0.0094,0.0754) θ2 0.0126 0.0104 0.0003 (0.0029,0.0416) PB β1 0.7498 0.1574 0.0039 (0.4968,1.1180) 688.9 GE β2 1.0510 0.2453 0.0058 (0.6632,1.6240) θ1 0.0063 0.0016 0.0000 (0.0036,0.0098) θ2 0.0080 0.0020 0.0000 (0.0045,0.0124) PB α1 0.0378 0.0247 0.0018 (0.0113,0.1058) 688.7 MW α2 0.0185 0.0170 0.0011 (0.0056,0.0695) β1 0.6436 0.1259 0.0096 (0.4083,0.9010) β2 0.7752 0.1457 0.0108 (0.4663,1.0110) φ1 0.0009 0.0007 0.0000 (0.0001,0.0026) φ2 0.0010 0.0008 0.0000 (0.0000,0.0031) copula functions reveal that, the estimates for each FGM bivariate distribution are similar to its corresponding product bivariate distribution. Furthermore, comparing between the fits of the bivariate distributions considering the FGM copula and the bivariate distributions considering the product copula reveal that, except for the FGMBMW distribution, the fits of the bivariate distributions considering the FGM copula are more efficient than the bivariate distributions considering the product copula function since they have the least information criteria values. Furthermore, the data is fitted to the FGMBNH distribu- 5 Aliyu & Usman / J. Nig. Soc. Phys. Sci. 5 (2023) 871 6 Table 3. posterior summaries considering FGM copula in the presence of co- variate - Kidney data model para- med- sd MC 95% DIC meter ian error CrI model λ 0.5223 0.3526 0.0145 (-0.3288,0.9705) 682.7 I φ10 -0.3115 0.2658 0.0111 (-0.8238,0.2320) φ11 0.5135 0.2207 0.0081 (0.0716,0.9167) φ12 -0.0015 0.0070 0.0002 (-0.0149,0.0126) φ20 -0.1204 0.2739 0.0107 (-0.6340,0.4389) φ21 0.0129 0.2110 0.0074 (-0.4173,0.4083) φ22 -0.0035 0.0077 0.0002 (-0.0181,0.0124) θ1 0.0137 0.0111 0.0004 (0.0035,0.0450) θ2 0.0111 0.0076 0.0002 (0.0033,0.0323) model λ 0.7428 0.2967 0.0052 (-0.1071,0.9901) 708.9 II φ1 0.2761 0.0417 0.0007 (0.2032,0.3684) φ2 0.3175 0.0541 0.0009 (0.2226,0.4331) θ10 -0.2323 0.3129 0.0046 (-0.8392,0.3709) θ11 0.2237 0.2903 0.0047 (-0.3455,0.7947) θ12 -0.0452 0.0111 0.0002 (-0.0662,-0.0226) θ20 -0.2378 0.3175 0.0053 (-0.8479,0.4009) θ21 0.0416 0.2903 0.0045 (-0.5374,0.6097) θ22 -0.0582 0.0105 0.0002 (-0.0780,-0.0371) model λ 0.7477 0.2645 0.0024 (0.0022,0.9893) 709.1 III φ10 -1.1140 0.1865 0.0039 (-1.4900,-0.7674) φ11 0.3020 0.2131 0.0022 (-0.1292,0.7034) φ12 0.0010 0.0080 0.0002 (-0.0138,0.0178) φ20 -0.9954 0.1906 0.0034 (-1.3930,-0.6387) φ21 0.0777 0.2234 0.0027 (-0.3676,0.5077) φ22 0.0040 0.0087 0.0002 (-0.0125,0.0218) θ10 -0.4131 0.3111 0.0038 (-1.0280,0.2000) θ11 0.1566 0.2959 0.0032 (-0.4234,0.7434) θ12 -0.0537 0.0135 0.0003 (-0.0803,-0.0273) θ20 -0.3647 0.3090 0.0030 (-0.9615,0.2536) θ21 -0.0040 0.2932 0.0030 (-0.5818,0.5714) θ22 -0.0687 0.0127 0.0002 (-0.0935,-0.0434) Table 4. posterior summaries considering FGM copula - Tobacco-stained- fingers data model para- med- sd MC 95% DIC meter ian error CrI FGM λ 0.9036 0.1220 0.0019 (0.5478,0.9967) 627.4 BNH φ1 0.2200 0.0723 0.0016 (0.1334,0.4149) φ2 0.2115 0.0446 0.0007 (0.1497,0.3258) θ1 1.4660 0.8177 0.0151 (0.4833,3.6510) θ2 1.4890 0.7025 0.0123 (0.6029,3.2950) FGM λ 0.9089 0.1155 0.0014 (0.5742,0.9966) 659.4 BE θ1 0.1417 0.0227 0.0003 (0.1017,0.1900) θ2 0.1037 0.0122 0.0001 (0.0817,0.1292) FGM β1 0.6686 0.0882 0.0009 (0.5051,0.8524) 631.1 BW β2 0.6192 0.0665 0.0009 (0.4964,0.7597) λ 0.9018 0.1244 0.0016 (0.5383,0.9963) θ1 0.2040 0.0369 0.0004 (0.1419,0.2863) θ2 0.2029 0.0325 0.0004 (0.1452,0.2723) FGM β1 0.6528 0.0994 0.0013 (0.486,0.8726) 634.2 BGE β2 0.5807 0.0742 0.0009 (0.4479,0.7384) λ 0.9032 0.1205 0.0015 (0.5487,0.9963) θ1 0.0733 0.0258 0.0003 (0.0333,0.1336) θ2 0.0504 0.0133 0.0002 (0.0282,0.0801) FGM α1 0.1898 0.0362 0.0006 (0.1273,0.2684) 634.4 BMW α2 0.1911 0.0325 0.0005 (0.1348,0.2607) β1 0.6302 0.0926 0.0014 (0.4597,0.8225) β2 0.5728 0.0745 0.0013 (0.4284,0.7207) φ1 0.0213 0.0255 0.0004 (0.0008,0.0962) φ2 0.0173 0.0177 0.0003 (0.0008,0.0662) λ 0.9010 0.1299 0.0022 (0.5132,0.9961) tion by taking sex and age as covariates. The covariates were Table 5. posterior summaries considering product copula - Tobacco-stained- fingers data model para- med- sd MC 95% DIC meter ian error CrI PB φ1 0.1917 0.0626 0.0010 (0.1173,0.3612) 643.4 NH φ2 0.2027 0.0443 0.0008 (0.1434,0.3138) θ1 1.6020 0.9106 0.0132 (0.5268,4.0580) θ2 1.5930 0.7913 0.0136 (0.6229,3.5410) PBE θ1 0.1277 0.0214 0.0002 (0.0907,0.1738) 678.7 θ2 0.1027 0.0123 0.0001 (0.0806,0.1288) PBW β1 0.6412 0.0873 0.0009 (0.4799,0.8255) 647.6 β2 0.6080 0.0675 0.0007 (0.486,0.7496) θ1 0.1896 0.0353 0.0004 (0.1287,0.2658) θ2 0.2030 0.0333 0.0004 (0.1452,0.2758) PB β1 0.6278 0.0967 0.0011 (0.4634,0.8398) 650 GE β2 0.5673 0.0747 0.0008 (0.438,0.7287) θ1 0.0587 0.0227 0.0002 (0.0244,0.1125) θ2 0.0476 0.0131 0.0001 (0.0263,0.0771) PB α1 0.1771 0.0351 0.0005 (0.1186,0.2570) 650.9 MW α2 0.1916 0.0329 0.0004 (0.1343,0.2625) β1 0.6025 0.0929 0.0013 (0.4298,0.7945) β2 0.5646 0.0752 0.0009 (0.4180,0.7148) φ1 0.0211 0.0249 0.0003 (0.0008,0.0941) φ2 0.0176 0.0180 0.0002 (0.0008,0.0673) assumed to have effect on φk, θk and on both φk and θk for k = 1, 2. The posterior summary statistics of the fits of the data are given in table 3. It is observed that the covariates have effect on φk. 3.2. Tobacco-stained-fingers data set In this section, we consider the tobacco stained finger data set obtained from [45]. The data consist of a sample of 143 smokers screened between March 2006 and January 2010 in a 180-bed community hospital in La Chauxde-Fonds, Switzer- land. Information on death and hospital admission of the pa- tients were collected until June 2014. For more details on this data set see [45]. To demonstrate the applicability of the introduced method- ology using this data set, it was considered as T1, the time before the first hospital readmission in smokers with stains on their fingers. On the other hand, the survival time of the patient with tobacco-tar stain on their fingers was considered as T2. This data is censored in case of death before the closure date. Similar procedure was followed in analyzing this data set as in section 3.1. The posterior summary statistics of the fits of the bivari- ate distributions to the tobacco-stained-fingers data are given in tables 4 and 5. The MC error of the posterior estimates of the fitted models considering the tobacco-stained-fingers data were all less than 120 of standard deviation of the estimates, this showed that the posterior estimates have reasonably good precision. Comparing the information criteria values in tables 4 and 5 showed that, the introduced BNH distribution fits the tobacco-stained-fingers data more efficiently than the compared bivariate distributions. Also, comparing between fits of the bivariate distributions considering the FGM copula and the bivariate distributions considering the product copula showed that, the bivariate distributions considering the FGM copula 6 Aliyu & Usman / J. Nig. Soc. Phys. Sci. 5 (2023) 871 7 fits the tobacco-stained-fingers data better than the bivariate distributions considering the product copula function. For the FGM bivariate models, the estimates of the dependent param- eter are similar for all the considered bivariate distributions. Moreover, the estimates of the shape and scale parameters of the bivariate distributions considering the FGM copula are similar to their bivariate distribution counterparts considering the product copula. 4. Conclusion In real life applications, lifetime data of two organs affected by a given type of treatment or therapy applied to the same in- dividual that present some type of weak dependence are usually obtained. In such a situation, the use of copula functions would be a very good alternative in modelling this type of bivariate lifetime data in presence of censored data and covariates. In this article, the analytical structures of the statistical methodol- ogy associated with the modelling of lifetimes assuming weak linear dependence was introduced. FGM copula function was used, however, there are many other copula functions that could be used to build new bivariate lifetime models depending on the assumptions established about the form of the relationship be- tween the bivariate lifetimes (T1 and T2). Also, marginal NH distributions were used since the distribution usually provide goodness of fit for lifetime data, as constant, decreasing or in- creasing hazard functions, features that are characteristic of the studied variables (T1 and T2). Hence, other univariate lifetime distributions such as inverted NH, unit NH, discrete NH distri- butions could also be used. Two real data sets: Kidney data and Tobacco-stained-fingers data sets were used in demonstrat- ing the applicability of the introduced methodology. The ob- tained results were compared assuming weak dependence that can be model using copula function, with those obtained con- sidering independence between lifetimes. Moreover, the intro- duced methodology is further compared with some competing bivariate distributions. The proposed approach could be useful in many areas of interest, including medical sciences, physical sciences and engineering studies. Acknowledgments We thank the referees for the positive enlightening com- ments and suggestions, which have greatly helped us in making improvements to this paper. References [1] M. H. Tahir, G. M. Cordeiro, S. Ali, S. Dey & A. Manzoor, “The inverted Nadarajah–Haghighi distribution: estimation methods and applications”, Journal of Statistical Computation and Simulation 88(14) (2018) 2775. [2] R. D. Gupta & D. Kundu, “Theory & methods: Generalized exponen- tial distributions”, Australian & New Zealand Journal of Statistics 41(2) (1999) 173. [3] S. Nadarajah & S. Kotz, “The beta exponential distribution”, Reliability engineering & system safety 91(6) (2006) 689. [4] S. Nadarajah & F. Haghighi, “An extension of the exponential distribu- tion”, Statistics, 45(6) (2011) 543. [5] A. Yakubu & S. I. Doguwa, “On the properties of the Weibull-Burr III distribution and its application to uncensored and censored survival data”, CBN Journal of Applied Statistics, 8(2) (2017) 91. [6] J. M. F. Carrasco, E. M. M. Ortega & G. M. Cordeiro, “A general- ized modified Weibull distribution for lifetime modeling”, Computational Statistics & Data Analysis 53(2) (2008) 450. [7] E. M. M. Ortega, G. M. Cordeiro & M. W. Kattan, “The log-beta Weibull regression model with application to predict recurrence of prostate can- cer”, Statistical Papers 54(1) (2013) 113. [8] M. Z. Raqab & D. Kundu, “Burr type X distribution: revisited”, Journal of probability and statistical sciences 4(2) (2006) 179. [9] A. I. Ishaq, A. Usman, M. Tasi’u, Y. Aliyu & F. A. Idris, “A new Weibull- Kumaraswamy distribution: Theory and applications”, Nigerian Journal of Scientific Research 16(2) (2017) 158. [10] I. Shah, B. Iqbal, A. M. Farhan, S. Ali & S. Dey, Unit Nadarajah and Haghighi distribution: properties and applications in quality control, Sci- entia Iranica, 2021. [11] S. Ali, S. Dey, M. H. Tahir & M. Mansoor, “The Poisson Nadarajah- Haghighi distribution: different methods of estimation”, Journal of Reli- ability and Statistical Studies (2021) 415 [12] M. Shafqat, S. Ali, I. Shah & S. Dey, “Univariate discrete Nadarajah and Haghighi distribution: Properties and different methods of estimation”, Statistica 80(3) (2020) 301-330 [13] U. Usman, S. Suleiman, B. M. Arkilla & Y. Aliyu, “Nadarajah-Haghighi model for survival data with long term survivors in the presence of right censored data. Pakistan Journal of Statistics and Operation Research 17(3) (2021) 695. [14] E. M. Almetwally, H. Z. Muhammed & E. S. A. El-Sherpieny, “Bivari- ate Weibull distribution: properties and different methods of estimation”, Annals of Data Science 7(1) (2020) 163. [15] X. Bai, Y. Shi, B. Liu, & Q. Fu, “Statistical inference of Marshall-Olkin bivariate Weibull distribution with three shocks based on progressive in- terval censored data”, Communications in Statistics-Simulation and Com- putation 48(3) (2019) 637. [16] E. A. El-Sherpieny, H. Z. Muhammed & E. M. Almetwally, “Fgm bi- variate weibull distribution”, In Proceedings of the Annual Conference in Statistics (53rd), Computer Science, and Operations Research, Institute of Statistical Studies and Research, Cairo University (2018) 55. [17] I. E. Gongsin & F. W. O. Saporu, “A bivariate conditional Weibull distri- bution with application”, Afrika Matematika 31 (2020) 565. [18] D. Kundu, & V. Nekoukhou, “On bivariate discrete Weibull distribution’, Communications in Statistics-Theory and Methods 48(14) (2019) 3464. [19] M. V. D. Peres, J. A. Achcar & E. Z. Martinez, “Bivariate modified Weibull distribution derived from Farlie-Gumbel-Morgenstern copula: a simulation study”, Electronic Journal of Applied Statistical Analysis, 11(2) (2018) 463. [20] M. El-Morshedy, M. S. Eliwa, A. El-Gohary & A. A. Khalil, “Bivariate exponentiated discrete Weibull distribution: statistical properties, estima- tion, simulation and applications”, Mathematical Sciences 14(1) (2020) 29. [21] J. A. Achcar, F. A. Moala, M. H. Tarumoto & L. F. Coladello, “A bivari- ate generalized exponential distribution derived from copula functions in the presence of censored data and covariates”, Pesquisa Operacional 35 (2015) 165. [22] R. Alotaibi, M. Khalifa, E. M. Almetwally & I. Ghosh, “Classical and Bayesian Inference of a Mixture of Bivariate Exponentiated Exponential Model”, Journal of Mathematics 2021 (2021) 1. [23] M. K. A. Elaal & R. S. Jarwan, “Inference of bivariate generalized ex- ponential distribution based on copula functions”, Applied Mathematical Sciences 11(24) (2017) 1155. [24] D. Kundu & R. D. Gupta, “Absolute continuous bivariate generalized exponential distribution”, AStA Advances in Statistical Analysis 95(2) (2011) 169. [25] S. M. Mirhosseini, M. Amini, D. Kundu & A. Dolati, “On a new abso- lutely continuous bivariate generalized exponential distribution”, Statisti- cal Methods and Applications 23 (2014) [26] S. M. Mirhosseini, M. Amini, D. Kundu & A. Dolati, “On a new abso- lutely continuous bivariate generalized exponential distribution”, Statisti- cal Methods & Applications 24(1) (2015) 61. [27] E. S. A. El-Sherpieny, E. M. Almetwally & H. Z. Muhammed, “Bayesian and non-bayesian estimation for the parameter of bivariate generalized 7 Aliyu & Usman / J. Nig. Soc. Phys. Sci. 5 (2023) 871 8 Rayleigh distribution based on clayton copula under progressive type-II censoring with random removal”, Sankhya A 2021 (2021) 1. [28] H. Z. Muhammed, “Bivariate inverse Weibull distribution”, Journal of Statistical Computation and Simulation 86(12) (2016) 2335. [29] D. Kundu & A. K. Gupta, “On bivariate inverse Weibull distribution”, Brazilian Journal of Probability and Statistics 31(2) (2017) 275. [30] A. S. Al-Moisheer, R. M. Alotaibi, G. A. Alomani & H. Rezk, “Bivari- ate mixture of inverse Weibull distribution: properties and estimation”, Mathematical Problems in Engineering 2020 (2020) 1. [31] M. S. Eliwa & M. El-Morshedy, “Bayesian and non-Bayesian estimation of four-parameter of bivariate discrete inverse Weibull distribution with applications to model failure times, football and biological data”, Filomat 34(8) (2020) 2511. [32] H. Z. Muhammed, & E. M. Almetwally, “Bayesian and non-Bayesian estimation for the bivariate inverse weibull distribution under progressive type-II censoring”, Annals of Data Science 10 (2023) 481. [33] S. Mondal & D. Kundu, “A bivariate inverse Weibull distribution and its application in complementary risks model”, Journal of Applied Statistics 47(6) (2020) 1084 [34] E. M. Almetwally & H. Z. Muhammed, “On a bivariate Frà c©chet distri- bution”, J Stat Appl Probab 9(1) (2020) 1. [35] H. Z. Muhammed, E. S. A. El-Sherpieny & E. M. Almetwally, ”Depen- dency measures for new bivariate models based on copula function”, In- formation Sciences Letters 10(3) (2021) 15. [36] S. Ali, M. Shafqat, I. Shah & S. Dey, “Bivariate discrete Nadarajah and Haghighi distribution: properties and different methods of estimation”, Filomat 33(17) (2019) 5589. [37] U. Usman & Y. Aliyu, “Bivariate Nadarajah-Haghighi distribution de- rived from copula functions: Bayesian estimation and applications”, Benin Journal of Statistics 5 (2022) 45. [38] D. Morgenstern, “Einfache beispiele zweidimensionaler verteilungen”, Mitteilingsblatt fur Mathematische Statistik, 8 (1956) 234. [39] D. J. G. Farlie, “The performance of some correlation coefficients for a general bivariate distribution”, Biometrika 47(3/4) (1960) 307. [40] E. J. Gumbel, “Bivariate exponential distributions”, Journal of the Amer- ican Statistical Association 55(292) (1960) 698. [41] P. K. Trivedi & D. M. Zimmer, ”Copula modeling: an introduction for practitioners”, Foundations and Trends R© in Econometrics 1(1) (2007) 1. [42] E. Z. Martinez, J. A. Achcar & T. R. Icuma, “Bivariate Basu-Dhar geo- metric model for survival data with a cure fraction”, Electronic Journal of Applied Statistical Analysis 11(2) (2018) 655. [43] D. J. Spiegelhalter, N. G. Best, B. P. Carlin & A. Van Der Linde, “Bayesian measures of model complexity and fit”, Journal of the royal statistical society: Series b (statistical methodology) 64(4) (2002) 583. [44] C. A. McGilchrist & C. W. Aisbett, ”Regression with frailty in survival analysis”, Biometrics, (1991) 461. [45] G. John, C. Louis, A. Berner & D. Genné, “Tobacco stained fingers and its association with death and hospital admission: A retrospective cohort study”, PloS one 10(9) (2015) e0138211. Appendices The following OpenBugs computational code can be used to obtained the posterior estimates from Bivariate Nadarajah- Haghighi distribution with parameters θ1, θ2, φ1, φ2 and λ. # The model assuming product copula function model { for (i in 1:N) { s1[i]< -exp(1-pow((1+theta1*t1[i]), phi1)) s2[i]< -exp(1-pow((1+theta2*t2[i]), phi2)) b1[i]<- pow((1+theta1*t1[i]), phi1-1) b2[i]<- pow((1+theta2*t2[i]), phi2-1) f1[i]< -theta1*phi1*b1[i]*s1[i] f2[i]< -theta2*phi2*b2[i]*s2[i] ft1t2[i]< - f1[i]*f2[i] del1[i]<- -f1[i]*s2[i] del2[i]<- -f2[i]*s1[i] st1t2[i]<- s1[i]*s2[i] L1[i]< - pow(ft1t2[i], d1[i]*d2[i])*pow(st1t2[i], (1-d1[i])*(1-d2[i]))*pow(del1[i], d1[i]*(1-d2[i]))* pow( del2[i], d2[i]*(1-d1[i])) L[i]<- abs(L1[i]) logL[i] < - log(L[i]) zeros[i] < - 0 zeros[i] ~ dloglik(logL[i]) } # prior distributions theta1 ~ dgamma(1,1) theta2 ~ dgamma(1,1) phi1 ~ dgamma(1,1) phi2 ~ dgamma(1,1) } # The model assuming FGM copula function model { for (i in 1:N) { s1[i]< -exp(1-pow((1+theta1*t1[i]), phi1)) s2[i]< -exp(1-pow((1+theta2*t2[i]), phi2)) b1[i]<- pow((1+theta1*t1[i]), phi1-1) b2[i]<- pow((1+theta2*t2[i]), phi2-1) f1[i]< -theta1*phi1*b1[i]*s1[i] f2[i]< -theta2*phi2*b2[i]*s2[i] f12[i]<- exp(2-pow((1+theta1*t1[i]), phi1)- pow((1+theta2*t2[i]), phi2)) ft1t2[i]< - theta1*theta2*phi1*phi2*f12[i]* b1[i]*b2[i]*(1+lambda*(1-2*s1[i])*(1-2*s2[i])) del1[i]<- f1[i]*s2[i]* (1+lambda*(1-2*s1[i]) *(1-s2[i])) del2[i]<- f2[i]*s1[i]* (1+lambda*(1-2*s2[i]) *(1-s1[i])) st1t2[i]<- f12[i]* (1+lambda*(1-s2[i])*(1-s1[i])) L1[i]< - pow(ft1t2[i], d1[i]*d2[i])*pow(st1t2[i], (1-d1[i])*(1-d2[i]))*pow(del1[i], d1[i]*(1-d2[i])) *pow( del2[i], d2[i]*(1-d1[i])) L[i]<- abs(L1[i]) logL[i] < - log(L1[i]) zeros[i] < - 0 zeros[i] ~ dloglik(logL[i]) } # prior distributions theta1 ~ dgamma(1,1) theta2 ~ dgamma(1,1) phi1 ~ dgamma(1,1) phi2 ~ dgamma(1,1) lambda~ dunif(-1,1) 8 Aliyu & Usman / J. Nig. Soc. Phys. Sci. 5 (2023) 871 9 } # The model assuming FGM copula function assuming covariate effect on the shape parameters model { for (i in 1:N) { phi1[i]<-exp(phi10+phi11*x1[i]+phi12*x2[i]) phi2[i]<-exp(phi20+phi21*x1[i]+phi22*x2[i]) s1[i]< -exp(1-pow((1+theta1*t1[i]), phi1[i])) s2[i]< -exp(1-pow((1+theta2*t2[i]), phi2[i])) b1[i]<- pow((1+theta1*t1[i]), phi1[i]-1) b2[i]<- pow((1+theta2*t2[i]), phi2[i]-1) f1[i]< -theta1*phi1[i]*b1[i]*s1[i] f2[i]< -theta2*phi2[i]*b2[i]*s2[i] f12[i]<- exp(2-pow((1+theta1*t1[i]), phi1[i])- pow((1+theta2*t2[i]), phi2[i])) ft1t2[i]< - theta1*theta2*phi1[i]*phi2[i]*f12[i]* b1[i]*b2[i]*(1+lambda*(1-2*s1[i])*(1-2*s2[i])) del1[i]<- f1[i]*s2[i]* (1+lambda*(1-2*s1[i]) *(1-s2[i])) del2[i]<- f2[i]*s1[i]* (1+lambda*(1-2*s2[i]) *(1-s1[i])) st1t2[i]<- f12[i]* (1+lambda*(1-s2[i])*(1-s1[i])) L1[i]< - pow(ft1t2[i], d1[i]*d2[i])*pow(st1t2[i], (1-d1[i])*(1-d2[i]))*pow(del1[i], d1[i]*(1-d2[i])) *pow( del2[i], d2[i]*(1-d1[i])) L[i]<- abs(L1[i]) logL[i] < - log(L1[i]) zeros[i] < - 0 zeros[i] ~ dloglik(logL[i]) } # prior distributions theta1 ~ dgamma(1,1) theta2 ~ dgamma(1,1) phi10~dnorm(0, 10) phi11~dnorm(0, 10) phi12~dnorm(0, 10) phi20~dnorm(0, 10) phi21~dnorm(0, 10) phi22~dnorm(0, 10) lambda ~ dunif(-1, 1) } In this codes, N is the sample size, S1[i] is the survival function given in equation (9), S2[i] is the survival function given in equation (10), f1[i] is the probability density function given in equation (7), f2[i] is the probability density function given in equation (8), ft1t2 [i] is the joint probability density function (pdf) given in equation (12), del1[i] is an expression given in equation (14), del2[i] is an expression given in equation (15), St1t2[i] is the joint survival function given in equation (11) and L[i] is the likelihood function given in equation (19). For the product copula function, the dependence parameter (lambda) assumes the value zero, and hence, the joint density function reduces to the product of the pdfs in equations (7) and (8). The joint survival function also reduces to the product of the survival functions in equations (9) and (10). For the models that assume covariate effect on the parameter(s) of the model, appropriate substitutions are made in the codes for the FGM copula function as discussed in the work. 9 Aliyu & Usman / J. Nig. Soc. Phys. Sci. 5 (2023) 871 10 Table 5. Kidney data Patient T1 T2 δ1 δ1 Sex Age Disease types 1 8 16 1 1 1 28 3 2 23 13 1 0 0 48 0 3 22 28 1 1 1 32 3 4 447 318 1 1 0 31.5 3 5 30 12 1 1 1 10 3 6 24 245 1 1 0 16.5 3 7 7 9 1 1 1 51 0 8 511 30 1 1 0 55.5 0 9 53 196 1 1 0 69 1 10 15 154 1 1 1 51.5 0 11 7 333 1 1 0 44 1 12 141 8 1 0 0 34 3 13 96 38 1 1 0 35 1 14 149 70 0 0 0 42 1 15 536 25 1 0 0 17 3 16 17 4 1 0 1 60 1 17 185 177 1 1 0 60 3 18 292 114 1 1 0 43.5 3 19 22 159 0 0 0 53 0 20 15 108 1 0 0 44 3 21 152 562 1 1 1 46.5 2 22 402 24 1 0 0 30 3 23 13 66 1 1 0 62.5 1 24 39 46 1 0 0 42.5 1 25 12 40 1 1 1 43 1 26 113 201 0 1 0 57.5 1 27 132 156 1 1 0 10 0 28 34 30 1 1 0 52 1 29 2 25 1 1 1 53 0 30 130 26 1 1 0 54 0 31 27 58 1 1 0 56 1 32 5 43 0 1 0 50.5 1 33 152 30 1 1 0 57 2 34 190 5 1 0 0 44.5 0 35 119 8 1 1 0 22 3 36 54 16 0 0 0 42 3 37 6 78 0 1 0 52 2 38 63 8 1 0 1 60 2 Table 5. Tobacco data n sex age T1 δ1 T2 δ2 control31 M 83 0.063014 1 0.013699 0 control34 M 36 6.69863 0 6.079452 0 control69 M 49 5.021918 0 4.40274 0 control55 M 39 5.59726 0 4.978082 0 control4 M 49 7.879452 0 3.876712 0 control39 F 85 5.309589 1 4.060274 0 control54 F 50 5.6 0 4.980822 0 control14 M 58 6.093151 1 2.490411 1 control24 F 80 0.065753 1 0.249315 0 control11 M 33 7.739726 0 0.646575 0 control5 M 64 1.405479 1 1.389041 1 control53 M 47 5.791781 0 3.4 0 control62 F 57 5.227397 0 4.545206 0 control27 M 46 7.435616 0 4.945206 0 control44 F 83 6.019178 0 1.654795 1 control17 M 64 0.668493 1 0.454795 1 control25 F 76 3.038356 1 0.227397 0 control16 M 64 0.347945 1 0.347945 0 control32 F 83 6.460274 0 2.753425 0 control57 M 71 5.536986 0 0.10137 1 control66 M 65 0.383562 1 0.375342 1 control8 M 50 7.605479 0 6.986301 0 control65 M 60 5.213699 0 4.594521 0 control59 F 41 5.536986 0 0.052055 0 control9 M 55 7.80274 0 7.183562 0 control33 F 35 6.506849 0 1.19726 0 control3 M 57 4.693151 1 2.139726 0 control63 M 55 0.575342 1 0.065753 0 control56 M 51 4.652055 1 4.109589 0 control7 F 51 7.80548 0 7.186301 0 control47 M 59 6.079452 0 2.767123 1 control71 M 50 4.978082 0 4.358904 0 control26 M 56 0.252055 1 0.252055 1 control15 M 70 0.005479 1 0.005479 0 control60 M 41 5.161644 0 4.542466 0 control36 F 82 2.736986 1 1.865753 0 control1 F 46 7.928767 0 0.654795 0 control68 M 74 5.410959 0 0.334247 1 control51 F 71 3.263014 1 3.126027 1 control46 F 46 0.257534 1 0.109589 0 control21 F 76 6.627397 1 2.161644 0 control41 F 59 6.243835 0 1.043836 1 control37 M 59 6.112329 0 3.271233 0 control18 F 58 7.835617 0 7.216438 0 control19 F 75 1.734247 1 0.69589 1 control67 F 68 5.367123 0 4.747945 0 control64 F 60 5.19726 0 4.578082 0 control45 F 71 6.172603 0 5.553425 0 control10 F 50 7.739726 0 7.120548 0 control61 M 64 5.254795 0 1.189041 0 control35 M 62 4.257534 1 2.553425 0 10 Aliyu & Usman / J. Nig. Soc. Phys. Sci. 5 (2023) 871 11 Table 5. Tobacco data cont. n sex age T1 δ1 T2 δ2 control13 M 37 3.512329 1 3.512329 0 control30 M 66 7.041096 0 1.252055 0 control70 F 83 0.167123 1 0.167123 0 control28 M 53 6.923288 0 0.052055 0 control2 M 50 7.969863 0 6.70411 0 control12 F 77 1.561644 1 1.49863 1 control58 F 81 5.534246 0 4.915069 0 control40 F 78 0.29589 1 0.167123 1 control29 M 65 6.926027 0 0.432877 1 control48 M 37 5.646575 0 0.638356 0 control43 M 64 0.243836 1 0.032877 1 control52 F 50 5.657534 0 5.038356 0 control49 M 76 0.682192 1 0.569863 0 control23 F 66 7.210959 0 0.216438 1 control20 M 66 7.827397 0 1.131507 0 control6 M 51 7.816438 0 0.276712 0 control50 M 40 5.619178 0 5 0 control22 F 65 7.128767 0 0.169863 0 control38 F 60 5.904109 0 3.413699 0 Tache40 M 68 1.183562 1 0.687671 1 Tache7 M 45 8.819178 0 1.747945 0 Tache64 M 54 7.419178 0 0.550685 0 Tache60 M 55 7.419178 0 1.59726 0 Tache61 F 71 7.512329 0 0.339726 0 Tache16 F 57 1.394521 1 0.071233 0 Tache57 M 64 0.315068 1 0.150685 1 Tache26 F 49 8.372602 0 0.463014 0 Tache73 M 65 2.180822 1 0.027397 0 Tache50 F 75 1.90411 1 1.69589 1 Tache42 M 73 7.260274 1 1.863014 1 Tache9 M 63 5.99726 1 0.482192 0 Tache25 M 62 0.128767 1 0.120548 1 Tache37 F 69 5.624658 1 0.320548 0 Tache6 M 50 8.819178 0 5.969863 0 Tache31 M 79 2.021918 1 0.076712 0 Tache67 F 88 1.780822 1 0.90411 0 Tache72 F 63 2.161644 1 0.071233 0 Tache8 M 71 8.873973 0 0.567123 1 Tache33 M 27 8.284931 0 7.665753 0 Tache30 M 56 8.265754 0 0.030137 0 Tache56 F 70 6.838356 1 2.616438 1 Tache13 M 69 7.882192 0 1.139726 1 Tache70 M 78 7.339726 0 1.936986 1 Tache19 M 50 0.679452 1 0.679452 0 Tache14 F 72 3.230137 1 0.8 0 Tache59 M 74 0.30137 1 0.30137 0 Tache22 M 52 2.123288 1 0.956164 1 Tache36 M 85 0.550685 1 0.139726 1 Tache53 M 50 3.60274 1 1.038356 0 Tache2 M 47 5.323287 1 4.482192 1 Table 5. Tobacco data cont. n sex age T1 δ1 T2 δ2 Tache18 M 69 8.452055 0 0.136986 0 Tache3 M 61 2.89589 1 1.649315 0 Tache1 M 67 0.071233 1 0.071233 0 Tache11 M 53 3.923288 1 1.994521 0 Tache35 F 44 8.30137 0 4.073973 0 Tache46 M 63 7.808219 0 0.09589 0 Tache47 F 66 7.769863 1 3.106849 1 Tache58 F 51 7.904109 0 0.424658 0 Tache21 F 69 2.232877 1 0.060274 0 Tache55 M 50 0.520548 1 0.109589 1 Tache29 M 45 0.052055 1 0.052055 0 Tache27 M 54 8.378082 0 1.032877 0 Tache5 M 38 8.821918 0 1.830137 0 Tache24 M 37 8.394521 0 7.775342 0 Tache51 M 62 7.547945 0 0.906849 0 Tache39 M 67 0.238356 1 0.013699 1 Tache54 M 53 7.616438 0 0.035616 0 Tache45 M 36 8.128767 0 7.509589 0 Tache20 M 49 8.380822 0 5.019178 0 Tache48 M 80 0.252055 1 0.016438 0 Tache4 F 64 0.136986 1 0.136986 0 Tache38 M 79 3.468493 1 0.39726 0 Tache49 M 67 7.654795 0 7.035616 0 Tache34 M 75 0.665753 1 0.386301 0 Tache65 F 63 7.128767 0 0.323288 0 Tache69 M 77 0.947945 1 0.071233 1 Tache17 M 83 6.306849 1 1.243836 0 Tache71 M 52 0.915069 1 0.868493 0 Tache41 M 59 7.846575 0 1.084931 0 Tache63 M 72 7.383562 0 0.213699 0 Tache15 M 50 1.20274 1 1.20274 0 Tache62 M 66 7.383562 0 5.243835 1 Tache12 F 75 0.021918 1 0.021918 0 Tache23 F 73 2.336986 1 2.336986 0 Tache66 M 44 7.189041 0 1.221918 0 Tache52 M 87 0.405479 1 0.405479 0 Tache10 F 69 2.386301 1 2.386301 0 Tache68 M 56 2.616438 1 2.610959 0 Tache28 M 54 8.40274 0 1.271233 1 Tache44 M 46 4.950685 1 1.991781 0 Tache43 F 69 8 0 1.30411 0 Tache32 F 71 5.134246 1 0.136986 1 11