Improved estimation of hazard function when failure information is missing not at random

ABSTRACT Hazard function plays a crucial role in survival analysis. Its estimation has garnered a lot of attention when the survival time variable suffers from right-censoring. Most of the existing works focus on the cases that failure information is complete or missing at random (MAR). When the censoring information is missing not at random (MNAR), statistical inferences on hazard function are very challenging. In this study, estimation of hazard function is addressed under the MNAR mechanism of the failure information. Three estimators are proposed by employing the techniques of correcting biases and adjusting weighting probabilities. These estimators are validated to be consistent and asymptotically normal. Simulation studies and two real-data analyses are performed to assess the proposed methods.


Introduction
Right-censored data arise in many research areas, such as clinical medicine, economics, reliability studies.Under the right-censored context, failure information is sometimes missing owing to many reasons (e.g.autopsies were not conducted or medical records were missing).Gijbels, Lin, and Ying (1993), McKeague and Subramanian (1998), van der Laan and McKeague (1998), and Gao and Tsiatis (2005) presented several typical examples of missing failure information.
The MAR mechanism means that the missingness is unrelated to the missing variable itself.The failure information is usually assumed to be missing at random, and a large body of literature focuses on estimation in survival models under this missing mechanism.More details have been reported in references Gao and Tsiatis (2005), Chen, He, Shen, and Sun (2009), Wang, Liu, and Liu (2009), Liu and Wang (2010), Song, Sun, Mu, and Dinse (2010), Lee, Cronin, Gail, Dignam, and Feuer (2011), Gouskova, Lin, and Fine (2017), Shen and Liang (2018), and relevant literature.Many researchers also focus on estimations in survival time models for homogeneous right-censored data.See Subramanian and Bean (2008a) and references therein.
Let T and C be survival and censoring time variables with distribution functions F and G, respectively.Assume that T and C are independent.One observes (X, δ, ξ) with X = min(T, C) and δ = 1(T ≤ C), where 1(A) is the indicator function of a set A, and ξ denotes the missing indicator with ξ = 1, when δ is observed and 0 otherwise.
For right-censored data, the missingness of the failure type is often related to the failure type itself.For example, the survival times for subjects withdrawn from the experiment or living at the end of the study (δ = 0) are known to be censored (ξ = 1).For a data set of 169 elderly women with stage II breast cancer, Sun, Xie, and Liang (2013) shows that the MAR assumption is rejected by conducting their testing procedure.
Clearly, by relaxing the independence assumption between the failure type and the missing indicator variable, the MNAR mechanism is more reliable.However, it is quite challenging to estimate the survival time functions under the non-ignorable mechanism since some key parameter is inestimable.To avoid the problem of inestimability, some additional information should be provided.
In this study, we aim to develop appropriate estimators of hazard function when failure type is missing not at random.To handle the censoring and the missing simultaneously, two parametric probability models are assumed.A conditional maximum likelihood method is employed to yield consistent estimators of the parameters in the probability models under the MNAR condition.Furthermore, estimators of the hazard function are built by combining the plug-in, inverse probability weighted, and kernel smoothing methods.Their consistency and weak convergence are then investigated rigorously.
The rest of this paper is organised as follows.In Section 2, a method is proposed to estimate the missing and censoring probability models.In Section 3, three estimators of the hazard function are built.Their asymptotic properties are studied in Section 4. In Section 5, an effective but simple bandwidth choice method is developed.Some simulation studies and two real-data analyses are conducted in Sections 6 and 7, respectively.The technical proofs of the main results are presented in the Appendix.Some preliminary lemmas and the detailed proofs are presented in the online Supplementary Material.

Bias-corrected imputation estimator
Imputation is very popular in missing data analyses.By imputing a reasonable value to replace the missing one, a completed data set is then constructed and then, an estimating procedure can be developed.Here, a completed data set {(X i , is a kernel function, and l n is a bandwidth sequence.
To estimate the hazard function, a natural imputation estimator is built as follows: where R i is the rank of X i , for i = 1, . . ., n, and k h n (x − X i ) is defined as h −1 n k((x − X i )/h n ) with a kernel function k(•) and a bandwidth h n .Under the MAR mechanism of the failure type, the above estimator tends to be consistent, as indicated by Wang, Dinse, and Liu (2012).However, under the MNAR mechanism, this imputation estimator is biased.The bias can be validated to be approximated by where f X (•) and H(•) are the density function and the cumulative distribution of X, respectively.This asymptotic bias vanishes when E(ξ |X = x, δ = 1) = E(ξ |X = x).Notably, the condition is exactly the MAR assumption.Furthermore, combining the prior information of the parametric model E(δ|X = x) = m(x, θ), a parametric imputation estimator is defined as follows: Its asymptotic bias under the MNAR mechanism is shown to be By correcting the bias of the estimator λp,I n (x), a bias-corrected imputation estimator can be constructed as where αn and θn are defined in Section 2, ûn , is an unbiased estimator of the asymptotic bias ABias( λp,I n (x)) in (3).

An adjusted complete-case estimator
Owing to a lack of experience in handling missing data, practitioners always ignore the missingness and apply a so-called complete-case estimator directly: , θ).Therefore, the complete-case estimator underestimates the true hazard function.
When the censoring information is fully observed, the estimator ) is consistent by summing over the observations with δ = 1.However, the complete-case estimator λC n (x) sums over only the observations with δ = 1 and ξ = 1.This causes a reduction of the summation terms and thus, underestimates the true hazard function.
By compensating the reduced summation terms on ξ = 0 in λC n (x), an adjusted complete-case estimator is defined as Note that {1 − π(X i , 1, αn )}m(X, θn ) is a consistent estimator of E((1 − ξ)δ|X).The introduction of this term just compensate the missing summation terms on ξ = 0 in λC n (x).

An adjusted inverse probability weighted (IPW) estimator
The IPW method is popular too in missing data analyses.Wang et al. (2012) proposed an IPW estimator of the hazard function: where ûn (x) and mn (x) are the Nadaraya-Watson kernel regression estimators of E(ξ |X = x) and E(δ|X = x) respectively, and have been defined in Section 3.1.The estimator λIP n (x) is shown to have desirable properties under MAR mechanism, such as double robustness.However, under the MNAR mechanism, the IPW estimator is asymptotically biased.The bias of the estimator λIP n (x) is asymptotically equal to Furthermore, considering the parametric model E(δ|X = x) = m(x, θ), a parametric IPW estimator can be defined as Under the MNAR mechanism, the estimator λp,IP n (x) is biased and its asymptotic bias can be validated to be Notice that under the MNAR mechanism, the following equation , an adjusted IPW estimator is proposed: The estimator λAI n (x) can also be derived by correcting the bias of λp,IP n (x).That is, Here, the second part of the right hand is an unbiased estimator of the asymptotic bias ABias( λp,IP n (x)) in (7).
Remark 3.1: As for the proposed three estimators, the problem of boundary effects should be considered.The method in Section 6.2 of Klein and Moeschberger (1997) can be employed.Specifically, on the interval [h n , τ − h n ], the estimators can be defined as shown above.Here, τ is a constant satisfying H(τ ) < 1.In the intervals [0, h n ] and [τ − h n , τ ], the kernel function should be modified.The scheme of modification can refer to monographs Klein and Moeschberger (1997) and Gasser and Müller (2006).Also, note that this is a common problem for the local kernel estimator of the hazard function.A lot of literature tends to ignore it.See Tanner and Wong (1983), Wang et al. (2012), and so on.

The asymptotic properties
In this section, the asymptotic properties of the above proposed estimators are established.In the subsequent parts, λ(x) denotes the true hazard function and λnj (x) represents λCI n (x) in (4), λAC n (x) in ( 6), and λAI n (x) in ( 8), for j = 1, 2, 3, respectively.Their uniform consistency and weak convergence are presented in Theorems 4.1 and 4.2, respectively.Theorem 4.1: Under Conditions (C1)-(C4) in Appendix 2, the following strong consistency holds: where τ is a constant satisfying H(τ ) < 1 and a.s.
The proofs of the above theorems are given in Appendix 3. Theorem 4.2 shows that the proposed three estimators obtain asymptotic variances as if the parameters θ and α were known.
From Theorem 4.2, the asymptotic mean integrated square errors (AMISE) can be calculated as Then the optimal bandwidths, which minimise the above AMISE j (h n ), are shown to be h oj = { σ 2 j (x) dx/ {λ (x)} 2 dx} 1/5 n −1/5 for the estimators λnj (x), j = 1, 2, 3.The bandwidth b n appears in the bias-corrected imputation estimator λCI n (x) and the adjusted IPW estimator λAI n (x).However, it doesn't affect AMISE j (h n ) in ( 9).The simulation results in Section 6 confirm that the bandwidth b n has little influence on λCI n (x) and λAI n (x).

The choice of bandwidths
Because all proposed estimators are local methods and only part of sample data are employed, the choices of bandwidths are critically important.The above optimal bandwidths are not feasible due to the unknown hazard function and its second derivative.In the following, a cross-validation least squares method is developed to generate practical bandwidths.
Consider the integrated squared error (ISE) functions , are distance measures between the estimators and the true function.Though the above unified representation of ISE is employed, the estimator λn2 (x) = λAC n (x) is only related to the bandwidth h n .Because the last term is unrelated to the bandwidth, the bandwidths are derived by minimising the following target functions Inspired by Patil (1993) and Wang et al. (2012), the bandwidths are chosen to minimise the computable objective function Here m is a large enough integer, and t 1 , t 2 , . . ., t m are grid points on [0, τ ].The superscript (−i) denotes the leave-one-out estimator.and 3, respectively.

Simulation studies
In this section, simulation studies are conducted to investigate the finite-sample properties of the proposed estimators.The following two settings are considered to generate observations.Setting 1.The survival time variable T and the censoring time variable C are generated from Weibull distributions with the hazard functions Setting 2. Let T and C generated from Gompertz distributions with the hazard functions In both settings, the missing indicator variable ξ is assumed to follow a probability model of the form π(x, δ, α) = 1/(1 + αx δ+1 ), α > 0. Clearly, the missing of the failure type variable δ is related to δ itself, which means that the observations are generated under the MNAR mechanism.Furthermore, the values of parameter θ in m(x, θ) and parameter α in π(x, δ, α) are specified to obtain certain censoring rate (CR) and missing rate (MR), which are shown in Table 1.
In the following, the mean integrated square error (MISE) of the proposed estimators will be calculated based on 1000 repetitions to evaluate the performance of the proposed estimators.For comparison, the MISE of the complete-case estimator λC n (x) in ( 5) is also computed.The full estimator, which is defined as by assuming that the failure type is fully observed, will act as a benchmark.We employ the biweight kernel function Estimation of the probability models: The parameters θ and α can be estimated by the conditional maximum likelihood method in Section 2. We calculate the bias (Bias), the sum of squares error (SSE), and the mean squared error (MSE) of estimators θn and αn under Settings 1 and 2 based on 1000 repetitions, and report the corresponding results in Table 2.We can observe that with the sample size n increasing, the values of the Bias, SSE, and MSE decrease, which verifies the consistency of the conditional maximum likelihood method.
Choice of bandwidths: Recall that the Bias-corrected imputation estimator λCI n (t) and the adjusted IPW estimator λAI n (t) depend on both h n and b n , whereas the adjusted completecase estimator λAC n (t) depends only on h n .Now we investigate the effects of h n and b n on the MISE of the proposed estimators.Specifically, MISEs are calculated for each pair (h n , b n ) in a 40 × 20 grid on the (0, 2] × (0, 1] plane.Figure 1  Furthermore, we investigate how h opt,n , the minimiser of the cross-validation function CV j (h n , b n ) in (10) for given b n , varies with b n .Here j = 1 and 3 correspond to λCI n (t) and λAI n (t), respectively.We choose 50 grids on (0, 1] for b n , and calculate the CVoptimal bandwidth h opt,n for each b n .Figure 3 displays the average values of h opt,n based In the follow-up simulations, we choose b n = 1.06 min{ σx , 0.75R x }n −1/5 by the rule of thumb (ROT), where σx is the sample standard deviation of X, and R x is the empirical inter-quantile range.The use of the empirical inter-quantile range aims for avoiding the influence of outliers.Then, the bandwidth h n is chosen by minimising (10) with the bandwidth b n replaced by the above ROT bandwidth.This operation can simplify calculations greatly.
Under the above-mentioned simulation configurations, the MISE of the proposed biascorrected imputation estimator (CI), the proposed adjusted complete-case estimator (AC), the proposed adjusted IPW estimator (AI), the complete-case estimator (CC), and the full estimator (FULL), are calculated based on 1000 repetitions.The results are reported in Table 3. From Table 3, it can be observed that all three proposed estimators yield relatively small MISEs than the complete-case estimator, and the MISE decreases with the increase of sample size or the decrease of CR.For the case that n = 200, CR = 20%, and MR = 20%, the proposed estimators λCI n (x) and λAC n (x) perform almost as well as the full estimator.Besides, the estimator λAI n (x) performs slightly worse than the other two estimators λCI n (x) and λAC n (x).In missing data analyses, the finite-sample property of the IPW-based estimator is generally a little worse than that of the imputation-based estimator.Note that the true hazard functions are λ(x) = 3x 2 and λ(x) = exp(x)/2 for Settings 1 and 2, respectively.We now plot the true hazard functions and average estimated curves of the five estimators for Settings 1 and 2 with sample size n = 60 in Figure 4. Figure 4 shows that the estimated curves of the three proposed estimators and the full estimator are very close to the true hazard curves, indicating that the proposed estimators perform well.However, the complete-case estimated curve deviates far from the true curve.
Performance under the MAR mechanism: We generate observations under Settings 1 and 2 except that the missing indicator variable ξ follows the parametric model π(x, α) = 1/(1 + αx), α > 0. The failure information δ is missing at random.The values of the parameter α vary in order to obtain the preset MR.We calculate the MISEs of the five  estimators and report the results in Table 4.It can be observed that the MISEs of the proposed estimators in Table 4 are comparable with those in Table 3, especially in the cases with larger sample sizes and lower CRs and MRs.It can be concluded that the proposed estimators are also feasible under the MAR mechanism.

Performance under misspecification:
A correctly-specified model for E(δ|X) = m(X, θ) is critical when we investigate the asymptotic properties of the proposed estimators.It is an interesting problem to explore the influence of model misspecification for E(δ|X) = m(X, θ).We generate observations under Setting 1 (Setting 2) but estimate the hazard functions by employing the parametric model m(X, θ) in Setting 2 (Setting 1).We rename the resultant setting as Setting 3 (Setting 4).Thus, the censoring probability models are incorrect.We plot the true hazard functions and estimated curves of the five estimators with sample size n = 60 in Figure 5 based on 1000 repetitions.From Figures 5 and 4, we can observe that the proposed estimators perform worse under misspecification.However, except for the adjusted IPW estimator λAI n (t) under Setting 4 with CR = 40%, the proposed estimators still perform better than the complete-case estimator that ignores the data missingness.

Analyses of primary biliary cirrhosis (PBC) data
In this subsection, we employ the proposed methods to analyse the PBC data.This data set was previously analysed by Subramanian and Bean (2008b) and is available at http://lib.stat.cmu.edu/datasets/pbc.The survival times of 418 patients with PBC are recorded in days.Among 418 patients, 144 patients died from the interesting disease PBC, and their survival times are observed accurately (ξ = δ = 1).Another 106 patients died from unknown causes, and their censoring indicators are missing (ξ = 0).The survival times of the left 168 patients are known to be censored and correspond to ξ = 1, δ = 0.  Here, we assume a censoring model m(x, θ) = θ 1 x/(1 + θ 2 x).If the survival time and the censoring time follow Weibull distributions, the above censoring probability model m(x, θ) holds.The response-probability model, π(x, 1, α), is assumed to be a kind of logistic model: α 1 exp(x)/{1 + α 2 exp(x)}.The kernel functions and the bandwidths are determined by the same methods in the simulation studies.The proposed three estimators and the complete-case estimator of the hazard function are calculated.The four estimated curves are plotted in Figure 6.
Figure 6 displays that the complete-case estimated curve is much lower than other estimated curves.This is consistent with the theoretical analysis that the complete-case estimator tends to underestimate the true curve.The estimated curves of the bias-corrected imputation and the adjusted complete-case estimators are close to each other.However, the adjusted IPW estimated curve deviates a little to the left from the bias-corrected imputation and the adjusted complete case estimated curves.
The estimated curves are all bell-shaped and reach their peaks between 3000 and 4000 days.It can be concluded that the PBC was progressively aggravated in ten years and then the risk of death owing to PBC decreases.Why the risk decreases is another problem worth exploring.Maybe, the PBC causes cancer ten years later, and then the patients would die from cancer rather than PBC.Anyway, the disease PBC needs to be treated as early as possible.

Real data analyses of non-renal vascular disease (NRVD) data
In this subsection, the proposed methods are applied to an NRVD data set analysed previously by Dinse (1986) and Wang et al. (2012).This data set reports the survival time and disease status at the death of 58 femaleFCFF mice from the results of the necropsy.Survival time was measured in days.Among 58 female mice, we considered 33 mice died with NRVD present.From the description of the data, it can be concluded that 8 mice died from NRVD(ξ = 1, δ = 1), 19 mice died from other known causes (ξ = 1, δ = 0), and 6 mice had an unknown cause of death (ξ = 0).We applied the proposed methods to estimate the hazard function for death due to NRVD.Analogously, the probability models m(x, θ) = θ 1 x/(1 + θ 2 x) and π(x, 1, α) = α 1 exp(x)/{1 + α 2 exp(x)} are employed.Other settings are the same as those in simulation studies.Figure 7 displays the estimated curves of the proposed three estimators and the complete-case estimator.
The curves are bell-shaped and reach peaks between 750 and 800 days.The curves of the proposed estimators are close to each other.However, the curve of the complete-case estimator is significantly lower than other curves.These results are consistent with the theoretical analyses.

Discussion
In this work, the problem of estimation of the hazard function is investigated when the censoring information is missing not at random.Three estimators are proposed and their asymptotic properties are established.The developed method dealing with MNAR failure type can be extended to estimate semi-parametric survival models.

Disclosure statement
No potential conflict of interest was reported by the author(s).Tanner, M.A., and Wong, W.H. (1983)

Appendices
Appendix 1. Notations , θ) , where ṁ and π denote the partial derivatives of m and π with respect to θ and α, respectively.
The conditional log-likelihood function in Equation (1) has a unique maximiser related to the parameter (θ, α).
Conditions (C1(i)) and (C1(ii)) are necessary for the parameter identifiability and the asymptotic expansion in Lemma A.1.Condition (C1(iii)) excludes the extreme cases that the three probability models are equal to 0 or 1. Condition (C2) is added to ensure that the conditional log-likelihood estimator is well-defined.Conditions (C3)-(C4) are common assumptions in nonparametric regressions.

Appendix 3. Proofs of Theorems
Four lemmas are first presented before proving the theorems.
with V presented in Appendix 1.
The above lemma can be proved by the routine method for a likelihood method and the details are omitted., δ, ξ) has the polynomial discrimination.The above lemma can be deduced from the theorem on Page 18 of Pollard (1984).

Lemma A.3: For the empirical distribution of X, H n (x), we have the following result:
The above lemma follows from the law of the iterated logarithm for an empirical distribution.
Lemma A.4: Under the same conditions of Theorem 4.2, for h n ≤ x ≤ τ − h n , the following result holds: Here, λ n (x) and σ 2 (x) are presented in Appendix 1.
The above lemma follows from Theorem 2 and Corollary 3 in Diehl and Stute (1988).

Proof of Theorem 4.1:
We only present the proof of the uniform convergence of the adjusted IPW estimator λAI n (x).Recalling the definition of λAI n (x) in ( 8), by some direct split calculations, it yields Equation ( 13) in Wang et al. (2012) indicates that sup From Lemma A.2, we have sup −→ 0. For any small constant ε > 0, when −→ 0. Similar to the proof of Theorem 2.6 in Li and Racine (2006), we can prove that sup From Lemma A.2, it is true that sup By Condition (C4) in Appendix 2, we obtain that sup The result that sup h n ≤x≤τ −h n |AI n24 (x)| a.s.
−→ 0 can be proved similarly.Therefore sup For AI n3 (x), we have The term AI n31 (x) can be decomposed into two parts: −→ 0, it can be proved that sup From ( A1)-(A5), Theorem 4.1 is then proved for λAI n (x).
Proof of Theorem 4.2: (I) Asymptotic expansion for λCI n (x).By some simple computations, we have By the facts that n − R i + 1 = n{1 − H n (X i −)} and Lemmas A.1-A.3, it is easy to validated that For the third term CI n3 (x), it can be split into four parts: By the results of Lemma A.1 and the law of large numbers, it is easy to prove that the last three main terms in (A8) are all o p (1).Therefore, it yields For CI n4 (x), we have s. and Lemma A.2, we obtain the following results: For CI n43 (x), we have By some tedious computations, we can validate that Therefore we can further prove that CI Then from (A10)-(A12), A9) and (A13), we have IF [1]  x (X i , δ i , ξ i ) + o p (1), (A14) with IF [1]  x (X, δ, ξ) presented in Appendix 1. (II) Asymptotic expansion for λAC n (x).Similar to (A6), we obtain that . For the third term AC n3 (x), we have The fourth main term in (A15) can be split into four parts: From (A15)-(A17), we have IF [1]  x (X i , δ i , ξ i ) + o p (1), (A18) with IF [1]  x (X, δ, ξ) presented in Appendix 1. (III) Asymptotic expansion for λAI n (x).By some split operations, we have =: Furthermore, the fourth term in (A19) can be decomposed into four parts: Similar to prove (A16), we can prove that AI n42 (x) = o p (1) and AI n43 (x) = o p (1).For the fourth term AI n44 (x), it can be decomposed into two parts: By the fact that sup h n ≤x≤τ −h n |1/ ûn (x) − 1/u(x)| = O( (nb n ) −1 log n + b 2 n ) a.s. and Lemma A.1, it is easy to obtain that the first term of AI n44 (x) is o p (1).Recalling that ûn (X i ) − u(X i ) = n j=1 {ξ j − u(X i )}γ , we have + o p (1) =: AI [1]  n44 (x) + AI [2]  n44 (x) + o p (1).Denote and then It can be proved that Further, by Condition (C1(i)) in Appendix 2 and the Taylor expansion, we can validate that AI [2]  n44 (x) = o p (1).So {ξ j − u(X j )}m(X j , θ) {1 − H(X j −)}u(X j ) k x − X j h n + o p (1), and then Thus we obtain that IF [2]  x (X i , δ i , ξ i ) + o p (1), (A20) with IF [2]  x (X, δ, ξ) presented in Appendix 1. (IV) The proof of the central limit theorem.From (A14), (A18), (A20), by the Lyapounov central limit theorem and Lemma A.4, Theorem 4.2 is then proved.

Figure 1 .
Figure 1.Simulation results of bandwidth selection.The MISE of λCI n (x) (upper row) and λAI n (x) (lower row) against bandwidths h n and b n for Setting 1 (left panel) and Setting 2 (right panel) with sample size n = 60, CR = 20%, and MR = 40%.

Figure 2 .
Figure 2. Simulation results of bandwidth selection.The MISE of λAC n (x) against bandwidth h n for Setting 1 (left panel) and Setting 2 (right panel) with sample size n = 60, CR = 20%, and MR = 40%.

Figure 3 .
Figure 3. Simulation results of bandwidth selection.The average CV-optimal bandwidth h opt,n of λCI n (x) (solid line) and λAI n (x) (dashed line) against bandwidth b n for Setting 1 (left panel) and Setting 2 (right panel) with sample size n = 60, CR = 20%, and MR = 40%.

Figure 5 .
Figure 5. Simulation results of hazard function estimation.Results under misspecification.True hazard function (thick line) and average estimated curves of the full estimator λF n (x) (thin line), λCI n (x) (dotted line), λAC n (x) (dash-dot line), λAI n (x) (thin dashed line), and λC n (x) (thick dashed line) for Setting 3 (left panel) and Setting 4 (right panel) with sample size n = 60.

Table 1 .
Values of parameters θ and α under different censoring rates and missing rates.
displays the MISE of λCI n (t) and λAI n (t) against h n and b n under Settings 1 and 2 with sample size n = 60, CR = 20%, and MR = 40%.Meanwhile, Figure 2 displays the MISE of λAC n (t) against h n with the same configuration.It can be seen clearly that b n is not critical for λCI n (t) and λAI n (t), whereas h n has a non-ignorable influence on all three estimators.

Table 2 .
Simulation results of estimators of parameters θ and α with different sample sizes, censoring rates, and missing rates.
on 1000 repetitions under Settings 1 and 2 with sample size n = 60, CR = 20%, and MR = 40%.We can observe that all the average CV-optimal bandwidths do not change much with b n , which indicates that b n are not critical for the minimisation of CV 1

Table 3 .
Simulation results of hazard function estimation.MISEs of the estimators with different sample sizes, censoring rates, and missing rates under the MNAR mechanism.

Table 4 .
Simulation results of hazard function estimation.
Notes: MISEs of the estimators with different sample sizes, censoring rates, and missing rates under the MAR mechanism.