Estimation on dependent right censoring scheme in an ordinary bivariate geometric distribution

ABSTRACT Discrete lifetime data are very common in engineering and medical researches. In many cases the lifetime is censored at a random or predetermined time and we do not know the complete survival time. There are many situations that the lifetime variable could be dependent on the time of censoring. In this paper we propose the dependent right censoring scheme in discrete setup when the lifetime and censoring variables have a bivariate geometric distribution. We obtain the maximum likelihood estimators of the unknown parameters with their risks in closed forms. The Bayes estimators as well as the constrained Bayes estimates of the unknown parameters under the squared error loss function are also obtained. We considered an extension to the case where covariates are present along with the data. Finally we provided a simulation study and an illustrative example with a real data.


Introduction
In the context of reliability and lifetime data analysis, lifetime variables are usually assumed to be nonnegative continuous random variables. Very often, lifetime values are available to the investigator in discrete form because they are measured discretely, that is, there is rounding to some small measurement unit. For instance, when we deal with the incubation period of a disease such as AIDS, or when the outcome variable of interest is the remission time of cancer, the time to event could be counted in number of weeks, months or so forth, see [24]. Sometimes the lifetime variables are discrete in theory meaning that time periods to the event of interest are countable instead of continuous. For instance, consider a study in which women who stop using oral contraception are followed until pregnancy. The outcome is defined as the number of menstrual cycles until pregnancy, see [15] for more examples.
In survival analysis, censoring is a key analytical problem that arises in a variety of ways. The most common type of censoring in both medical and engineering experiments is the right censoring. Statistical inference based on right censorship has been discussed extensively in survival analysis literature. For more details see, [1,4,16,22] and references therein.
Despite the high frequency of discrete measuring of lifetime values, in both clinical trials and engineering experiments, just a few works addressed the censoring in discrete setup. Rezaei and Arghami [31], obtained the maximum likelihood estimators (MLEs) of the parameters while the lifetime variables has an Exponential-geometric distribution in the presence of right censoring, Davarzani and Parsian [10], proposed the Middle censoring in a discrete setup and Davarzani and Parsian [11], considered an independent right censoring when both the lifetime and censoring variables have geometric distributions. Jammalamadaka and Leong [17] considered analysis of discrete lifetime data under middle-censoring in the presence of covariates.
Usually, lifetime variables are assumed to be independent from the censoring variables. Whereas, there are many situations that the lifetime variable could be dependent on the time of censoring. For instance, in some clinical trials, potentially harmful therapies (e.g. radiation or chemotherapy) may have side effects on patients, depending on the toleration of patients and the harmful degree of the therapy. Then patients may need to be withdrawn from the clinic. In such cases the time of withdrawal depends on the risk of time-to-event (e.g. death). In some studies, clinicians are interested to see the time to death from a particular disease (main disease) but the patients died from another competing disease that could be correlated with the main disease. In this case the event of interest (time to death from the main disease) is censored at the time to death from the competing disease. In industrial experiments, it may occur that a component is taken away from the experiment (is censored) because it shows some sign of future failure.
Recently, dependent censoring has been discussed in a few works in the literature in continuous setup. In most of them the dependence between the censoring and lifetimes variable modeled by allocating a known bivariate distribution to the lifetime and censoring variables or by applying a Copula functions. Jiang et al. [18] considered the problem of right censoring in a semi-parametric model where the dependence between the censoring mechanism and lifetimes is modeled by a gamma frailty copula. Braekers and Veraverbeke [6], obtained a copula-graphic estimator for the conditional survival function. Davarzani and Parsian [9], considered dependent right censoring where the lifetime and censoring variables have a Marshall-Olkin bivariate exponential (MOBE) distribution. Davarzani et al. [12], dealt with a dependent middle censoring in a dependent setup using MOBE distribution.
In this paper we propose the dependent right censoring scheme in discrete setup using a discrete bivariate distribution. Several bivariate discrete models have been proposed in statistics to represent lifetime data, which can be applied in this regard. Basu and Dhar [5] proposed a bivariate geometric distribution. Arnold [2] introduced multivariate geometric distributions which showed that it leads in a natural way to the Marshall-Olkin multivariate exponential distribution [28]. Nair and Nair [29] studied the characterizations of the geometric distributions. Basu and Dhar [5] introduced a bivariate geometric model which is analog to the bivariate distribution of Marshall and Olkin [27]. In our proposed scheme, it is not possible that a response of an individual to be observed and censored at the same time. Hence, among all the proposed geometric distributions, we found the bivariate geometric distribution proposed by Arnold [2] more applicable and flexible in the presence of dependent censoring of such an scheme.
The rest of the paper is organized as: In Section 2, we briefly introduce the bivariate geometric distribution provided by Arnold [2]. In Section 3, we introduce the potential data under the right censoring scheme and derive the likelihood and obtain the MLEs of the unknown parameters with their risks in closed forms. In Section 4, we deal with Bayes estimates of the unknown parameters and survival function under the squared error loss (SEL) function. Using constrained Bayes methodology, we obtained the constrained Bayes estimates of the unknown parameters under the SEL function. In Section 5, we propose the Bayes prediction of the next observation based on the observed data. In Section 6, we consider estimation under dependent right censoring in the presence of covariates. We obtain MLEs of the unknown parameters via method of scoring. Finally in Section 7, we provide a simulation study and an illustrative example in connection with our obtained estimators using a simulation data and a real data.

Preliminaries
Consider a sequence {T 1 , T 2 , T 3 , . . .} of i.i.d. random variables with states S 1 , S 2 and S 3 , that refer to event of type one, type two and type three, respectively. In the context of survival analysis S 1 , S 2 and S 3 states, can be supposed as failure (as an event of interest), censoring and neither censoring nor failure, respectively. Define the corresponding probability of every state as follows where, θ 1 , θ 2 , θ 3 > 0 and θ 3 = 1 − θ 1 − θ 2 . Define X and Y as X = min{x : T x ∈S 1 }, Y = min{y : T y ∈S 2 }, then we obtain the joint probability mass function (pmf) of X and Y as follows: see Omey and Minkova [30] for more details. We call this distribution 'Ordinary Bivariate Geometric' (OBG) with parameter θ = (θ 1 , θ 2 ) T , and denote it by OBG(θ 1 , θ 2 ). The marginal distributions of X and Y are ordinary univariate geometric distributions with parameters θ 1 and θ 2 , respectively.
Proof: If V = 0, the joint pmf of Z and V is given by If V = 1, the pmf is given by Now, combining Equations (2) and (3), we obtain the joint pmf of (Z, V) as follows: Now it is easy to verify that the pmf of (Z, V) in Equation (4) is a multiplication of marginal pmfs of Z and V , and this completes the proof.

MLEs of θ 1 and θ 2 in the presence of right censoring
Suppose that we have a random sample of n individuals and we only check them once in a season. Although, due to problems such as lack of follow up and competing risks, we do not always have opportunity of observing X 1 , . . . , X n . The corresponding potential censoring times are denoted by Y 1 , . . . , Y n . Hence the observed data are . . , n and has pmf obtained in Equation (4). This kind of censoring is called right censoring, see Davarzani and Parsian [9]. In this setup we assume the ith patient has lifetime variable X i that depends on the right censoring variable, Y i and (X i , Y i ), i = 1, . . . , n, have an OBG(θ 1 , θ 2 ) distribution defined in Equation (1). Now using Equation (4), the likelihood function based on actual observed data w is given by where Thus the MLEs of θ 1 and θ 2 , are the solutions of the following log-likelihood equations: And it is also easy to verify that the survival function of X i is given by Notice that from Lemma 2.1, we know that Z ∼ NB(n, θ 1 + θ 2 ) and V ∼ Bin(n, θ 1 / (θ 1 + θ 2 )) are independent. It is straightforward to show that V/n converges in distribution to a Normal distribution and Z/n converges in probability to 1/(θ 1 + θ 2 ). Following the well known asymptotic properties of MLEs, we calculated the Fisher information of θ and we discovered that √ n(θ − θ) converges in distribution to Bivariate Normal with mean 0 and variance covariance matrix The following Lemma will help us to obtain the risk functions of MLEs in a closed form. Lemma 3.1: Let Z has a NB(n, θ 1 + θ 2 ) distribution, then we have: Proof: See the Appendix (see supplementary material). Now, from Lemma 3.1, and using the facts that E(V) = nθ 1 /(θ 1 + θ 2 ) and E(V 2 ) = nθ 1 (θ 2 + nθ 1 )/(θ 1 + θ 2 ) 2 , the risk functions ofθ 1 andθ 2 under the SEL function, can be obtained as follows:

Bayesian analysis
It is well known, see [23], that estimating with uniformly minimum risk do not exist and there are two different ways in obtaining optimal estimates. One way is to restrict the estimators by requiring to satisfy some condition enforcing a certain degree of impartiality. Unbiasedness and invariance are two different impartiality conditions of this method. However, this approach of minimizing the risk uniformly in θ after restricting the estimators to be considered has limited applicability. Alternative and more general approach is to minimize an overall measure of the risk function associated with an estimator without restricting the estimators to be considered. One global measure of the size of the risk associated with an estimator of θ leads to Bayes estimators. Although we derived the MLEs of the parameters and adopting an asymptotic theory to describe the variation of the derived estimations, in the censored data setting, one problem could be the scarcity of data. Therefore, the Bayesian approach provides a coherent framework. In this section, we consider the Bayes estimation and constrained Bayes estimation of unknown parameters of interest in Sections 4.1 and 4.2, respectively.

Bayes estimation
We believe, as stated in [3], that from a strictly Bayesian viewpoint there is clearly no way in which one can say that one prior is better than any other. Presumably one has one's own subjective prior and must live with all of its lumps and bumps. It is more frequently the case that we elect to restrict attention to a given flexible family of priors and we choose one from that family that seems to best match our personal beliefs. If all members of a family of prior distribution possess undesirable features, such as insensitivity to the data, we should elect to use a different family of priors. Krishna and Pundir [21] obtained the Bayes estimate of parameters for some bivariate geometric distribution using Dirichlet distribution as prior. Li and Dhar [24] used a uniform prior distribution instead of using Dirichlet distribution.

Constrained Bayes estimation
It is well known that the Bayes estimation of θ under the SEL function is the vector of posterior means. As pointed out in [13,14,25], the empirical histogram of the posterior means is underdispersed as an estimate of the histogram of the corresponding population parameters. Thus adjustment of Bayes estimators is needed in order to meet the twin objective of simultaneous estimation and closeness of the histogram of the estimates with the posterior estimates of the parameter histogram. Hence, we consider constrained Bayes estimations of θ 1 and θ 2 under the SEL function by matching the first two empirical moments of the set of estimates with the corresponding posterior moments of the population parameters.
In this subsection, we will derive the constrained Bayes (CB) estimateθ within the class of all estimates t ≡ t(w) = (t 1 (w), t 2 (w)) of θ that satisfy Following [13], theθ CB is given bŷ | w), where 1 2 is a 2component column vector with each element equal to one, I 2 is the identity matrix of order 2 and J 2 = 1 2 1 T 2 . When θ has a Dir 2 (α 0 , α 1 , α 2 )-prior distribution, we have and H 1 (w) and H 2 (w) reduce to

Bayes prediction
In this section we obtain the Bayes prediction of W n+1 based on the observed data W = (W 1 , . . . , W n ). The predictive distribution of W n+1 | W = w is given by We are interested to predict Z n+1 under the SEL function and V n+1 under zero-one loss function. These loss functions are chosen based on the nature of the variables that must be predicted. It is obvious that E(Z n+1 | W = w) is the Bayes prediction of Z n+1 and reduces toZ where n > 2.
On the other hand the Bayes prediction of V n+1 under zero-one loss function is given by the posterior mode which is called the maximum a posteriori (MAP) prediction. For this purpose the marginal predictive distribution of V n+1 | W = w is It is easy to verify that the MAP predictor of V n+1 is given bỹ

The model in the presence of covariates
In this section, we consider a geometric lifetime with dependent right censoring in the presence of covariates. Let U T i = (U i1 , U i2 , . . . , U im ) be the covariate vector assigned to an individual with lifetime X i , that is, each person has a survival time X i and covariates U i . We assume ln{θ 1i /(1 − θ 1i )} = β T U i , where θ 1i denotes the probability of failure for individual i, i = 1, 2, . . . , n. Further, assume β = (β 0 , β 1 , . . . , β m ) is an unknown vector of coefficients. To set an appropriate model with substituting θ 1i = (exp{β T U i })/(1 + exp{β T U i }) in Equation (5), the likelihood function of β and θ 2 is obtained as follows: Hence, the log-likelihood function is and the MLEs of β and θ 2 can be obtained from the following equations: From the above equations we can't find the MLEs of the parameters in closed forms. However, they can be solved numerically. The second derivatives of log-likelihood are given by 2 , For convention let β * T = (β T , θ 2 ). Therefore, the Hessian of the log-likelihood matrix has the following form: Now using the method of scoring, at (j + 1)th stage we reach at Keep repeating this procedure till some convergence criterion for terminating the iteration obtain, that is, repeat the process until differences between successive estimates are sufficiently close to zero. The resulting estimates are denoted byβ * = (β,θ 2 ). Under regularity conditions,β * converges in distribution to multivariate normal distribution with mean equal to true parameter value and variance-covariance given by the inverse of the information matrix, that is,β * has approximately N m+1 (β * , I −1 (β * )).

Numerical exploration
In this section we conduct a simulation study and analyze an illustrative example to demonstrate the application of the proposed approach.

An illustrative example
In this section, we apply the proposed estimators and regression model on a data set about the survival times of patients with prostate cancer in stages 3 and 4 that have been followed-up during 1967-1969, see [7]. This data set has been analyzed in several studies, for example, [8,19,26]. In the original data set, patients could die from prostate cancer, cardiovascular disease (CVD), or other causes. In this study, we just include 161 patients (from 502 patients), who were died from prostate cancer or CVD. We are interested to analyze the time of death only from prostate cancer in the absence and in the presence of some covariates. These covariates are labeled by RX, AG, WT, PF, HX, HG, SZ, and SG, correspond, to treatment group (0, receiving low-dose of diethylstilbestrol (DES) and 1, receiving high-dose of DES), age group, weight index, performance rating, history of Table 2. Estimates (SREMSE) of θ 1 = 1/2 and θ 2 = 1/3 with α 0 = 1/3, α 1 = 1 and α 2 = 2/3.   CVD, serum hemoglobin, size of primary lesion, and Gleason stage/grade category. We had one missing data in the SZ column and four missing data in the SG column, which we imputed them to the nearest integer of their column mean, respectively. In this study, time to death from prostate cancer has been taken as the lifetime of interest that can be censored at time to death from CVD (as a competing risk). In this case we can never fully determine whether time to death from prostate cancer and time to death from CVD (as the censoring cause) are independent, see chapter 9 of [20]. Using the above information, it is easy to find the values of z = 4306 and v = 130, that depicts 80.75% of patients died from prostate cancer and 19.25% died from CVD and censored. Let (α 0 , α 1 , α 2 ) = (1/6, 1/2, 1/3), then the MLE, Bayes and Constrained Bayes Estimates of (θ 1 , θ 2 ) are (0.03019043, 0.007199257), (0.03029951, 0.007274979) and (0.03034977, 0.007224719), respectively. Figure 1 displays the ML estimated survival curves for patients in different groups of treatment that indicates, the estimated survival curve for the group with low-dose of DES lies below that for the group with high-dose of DES. For testing differences in survival curves, that is, H 0 : S 1 (t) = S 2 (t) vs H 1 : S 1 (t) = S 2 (t) we split the data into two groups namely, the group with low-dose of DES with n L patients and parameters θ 1L , θ 2L and the group with high-dose of DES with n H patients and parameters θ 1H , θ 2H . Now using the asymptotic distributions of the ML estimates of θ 1 and θ 2 in these two group, that ) and x is the integer part of x. In Byar sub-data, T = 1.236394, c 2 = 0.5793856, df = 79 and approximate p−value = 0.215 which guarantees acceptance of the null hypothesis H 0 : S 1 (t) = S 2 (t) at 10% significance level.
To analyze the effect of covariates in the Byar sub-data, we entered eight covariates to the model. The outcome of running the data with eight covariates show that estimates of the regression coefficients are very close to zero except for RX, AG, WT, HX and HG with estimated coefficients (standard errors of coefficients) −0.30786745(0.188251745), −0.01311510(0.007513652), −0.00866909(0.006332463), −0.29698126(0.190059084) and −0.110739786(0.042973908). Calculating |β/s.e.(β)| leads to 1.635403, 1.745503, 1.368992, 1.562573 and 2.576908 that guarantees the significance of AG at 10% significance level and HG at the 5% significance level. Kay [19] at his model (2) dealt with the effects of the explanatory variables on the causes separately. He investigated that RX, HG, SZ and SG affect risk of cancer death with higher values being associated with higher risk. Since our subpopulation is different from Kay [19], we discovered that at the 5% significance level HG is significant.

Discussion
Measuring lifetime values discretely or rounding them to small measurement units is very common clinical investigations. Depending on the follow up period of each individual, these values may be either real lifetimes or censored lifetimes. Usually, lifetime variables are assumed to be independent from the censoring variables, while there are many situations that the lifetime variable could be dependent on the time of censoring. In this paper we considered dependent right censoring while the lifetimes are in a discrete form. We proposed Maximum Likelihood, Bayes and Constrained Bayes estimators of the survival function and unknown parameters in closed forms. The simulation study demonstrated that all three methods perform well, even for small sample sizes, with respect to the SREM-SEs and biases of estimations. However, the Bayes estimators outperform with respect to the SREMSE.
A multivariate regression model provided in this setup and applied to analyzing a real data. To the best knowledge of the authors there exists no key regression model dealing with dependent discrete censored data for assessing the validity of the results.
Authors observed that there are some discrepancies between the results obtained from the proposed model and results obtained from Cox regression model in terms of obtained p-values; however, they are not (very) inconsistent especially in terms of magnitudes and signs of estimated parameters. This inconsistency might be a consequence of ignoring the correlation between the lifetime and censoring variables in using Cox regression models.
Although, it will be difficult, but it might be interesting to consider a bivariate discrete distribution with nonconstant hazard rate which might be more fitted to clinical data. We believe, more work is needed along these directions.

Disclosure statement
No potential conflict of interest was reported by the authors.

Funding
Ahmad Parsian's research was supported by the University of Tehran.