# Sampling with censored data: a practical guide

Pedro L. Ramosa, Daniel C. F. Guzmana,b, Alex L. Motaa,b,
Francisco A. Rodriguesa and Francisco Louzadaa

aInstitute of Mathematics and Computer Science, University of São Paulo, São
Carlos, Brazil
bDepartment of Statistics, Federal University of São Carlos, São Carlos,
Brazil Corresponding author. Email<EMAIL_ADDRESS>
###### Abstract

In this review, we present a simple guide for researchers to obtain pseudo-
random samples with censored data. We focus our attention on the most common
types of censored data, such as type I, type II, and random censoring. We
discussed the necessary steps to sample pseudo-random values from long-term
survival models where an additional cure fraction is informed. For
illustrative purposes, these techniques are applied in the Weibull
distribution. The algorithms and codes in R are presented, enabling the
reproducibility of our study. Key Words: Simulation; Censored data; Sampling;
Cure fraction, Weibull distribution.

##  Introduction

The recent advances in data collection have provided a unique opportunity for
scientific discoveries in medicine, engineering, physics, and biology (Shilo
et al., 2020; Mehta et al., 2019). Although some applications offer a vast
amount of data, as in astronomy or genetics, there are some other scientific
areas where missing recordings are very common (Tan et al., 2016). In many
situations, the data collection is limited due to certain constraints, such as
temporal or financial restrictions, which do not allow the acquisition of the
data set in its complete form. Mainly, missing values is very common in
clinical analysis, where, for instance, some patients do not recover from a
given disease during the experiment, representing the missing observations.

Survival analysis, which concerns the modeling of time to event data, provides
tools to handle missing data in experiments (Klein and Moeschberger, 2006).
Different approaches to deal with partial information (or censored data) have
been discussed through the survival analysis theory, where the results are
mostly motivated by time-to-event outcomes, widely observed in applications,
including studies in medicine and engineering (Klein and Moeschberger, 2006).
Although most books and articles have focused on discussing the theory and
inferential issues to deal with censored data (Bogaerts et al., 2017), it is
challenging to find comprehensive studies discussing how to obtain pseudo-
random samples under different types of censoring, such as type I, type II,
and random censoring. Type II censoring is usually observed in reliability
analysis where components are tested, and the experiment stops after a
predetermined number of failures (Singh et al., 2005). Hence, the remaining
components are considered as censored (Tse et al., 2000; Balakrishnan et al.,
2007; Kim et al., 2011; Guzmán et al., 2020). On the other hand, type I
censoring is widely applied in (but not limited to) clinical experiments where
the study is ended after a fixed time, and the patients who did not experience
the event of interest are considered as censored (Algarni et al., 2020).
Finally, the random censoring consists of studies where the subjects can be
censored during any experiment period, with different times of censoring (Qin
and Jing, 2001; Wang and Wang, 2001; Wang and Li, 2002; Li, 2003; Ghosh and
Basu, 2019). Random censoring is the most commonly used in scientific
investigations since it is a generalization of types I and II censorings.

This paper provides a review of the most common approaches to simulate pseudo-
random samples with censored data. A comprehensive discussion about the
inverse transformation method, procedures to simulate mixture models, and the
Metropolis-Hastings (MH) algorithm is presented. Further, we provide the
algorithms to introduce artificially censoring, including types I, II, and
random censoring, in pseudo-random samples. The approach is illustrated in the
Weibull distribution (Weibull, 1939), and the codes are present in language R,
such that the interested readers can reproduce our analysis in other
distributions. Additionally, we also present an algorithm that can be used to
generate pseudo-random samples with random censoring and the presence of long-
term survival (see (Rodrigues et al., 2009) and the references therein). This
algorithm will be useful to model medical studies where a part of the
population is not susceptible to the event of interest (Chen et al., 1999;
Peng and Dear, 2000; Chen et al., 2002; Yin and Ibrahim, 2005; Kim and Jhun,
2008).

This paper is organized as follows. Section $2$ discuss the main properties of
the Weibull distribution, which is the most well-known lifetime distribution.
Section $3$ describes the most common approaches to generate pseudo-random
samples. Section $4$ presents a step-by-step guide to simulate censored data
under different censoring mechanisms and cure fractions. The parameter
inference to estimate the parameters assuming different schemes to validate
the proposed approaches is discussed in section $5$. Section $6$ presents an
extensive simulation study using different metrics to validate the results
under many samples to illustrate the censoring effect in the bias and the mean
square error of the estimator. Finally, section 7 summarizes the study and
discusses some possible further investigations.

##  Weibull distribution

Initially, we describe the Weibull distribution, that will be considered in
the derivation of proposed sampling schemes. Let $X$ be a positive random
variable with a Weibull distribution (Weibull, 1939), then its probability
density function (PDF) is given by:

$f(x;\beta,\alpha)=\alpha\beta x^{\alpha-1}\exp\left\\{-\beta
x^{\alpha}\right\\},\quad x>0,$ (1)

where $\alpha>0$ and $\beta>0$ are the shape and scale parameters,
respectively. Hereafter, we will use $X\sim\text{Weibull}(\beta,\alpha)$ to
denote that the random variable $X$ follows this distribution. The
corresponding survival and hazard rate functions of the Weibull distribution
are given, respectively, by

$\displaystyle S(x;\beta,\alpha)=\exp\left\\{-\beta
x^{\alpha}\right\\},\quad\mbox{and}\quad
h(x;\beta,\alpha)=\dfrac{f(x;\beta,\alpha)}{S(x;\beta,\alpha)}=\alpha\beta
x^{\alpha-1}.$ (2)

The survival function of the Weibull distribution is proper, that is, $S(0)=1$
and $S(\infty)=\lim_{x\to\infty}S(x)=0$. Additionally, its hazard rate is
always monotone with constant ($\alpha=1$), increasing ($\alpha>1$) or
decreasing ($\alpha<1$) shapes.

The Weibull distribution has two parameters, which usually are unknown when
dealing with a sample. The maximum likelihood approach is the most common
method used to estimate the unknown parameters of the model. In this case, we
have to construct the likelihood function defined as the probability of
observing the vector of observations $\bm{x}$ based on the parameter vector
$\bm{\theta}$. We obtain the maximum likelihood estimator (MLE) by maximizing
the likelihood function. The likelihood function for $\bm{\theta}$ can be
written as

$L(\bm{\theta;x,\delta})=\prod_{i=1}^{n}f(x_{i}|\bm{\theta})^{\delta_{i}}S(x_{i}|\bm{\theta})^{1-\delta_{i}},$

where $\delta_{i}$ is the indicator of censoring in the samples with
$\delta=1$ if the data is not censored and $\delta=0$ otherwise. Note that
when $\delta_{i}=1,i=1,\ldots,n$ the likelihood function reduces to
$L(\bm{\theta,x})=\prod_{i=1}^{n}f(x_{i}|\bm{\theta})$. The evaluation of the
partial derivatives of the equation above is not an easy task. To simplify the
calculations, the natural logarithm is considered because the maximum value at
log likelihood occurs at the same point as the original likelihood function.
Notice that the natural logarithm is a monotonically increasing function. We
consider the $\log$-likelihood function that is defined as
$\ell(\bm{x,\delta;\theta})=\log L(\bm{\theta;x,\delta})$. Furthermore, the
MLE, defined by $\hat{\bm{\theta}}(\bm{x})$, is such that, for every $x$,
$\widehat{\bm{\theta}}(\bm{x})=\arg\max_{\bm{\theta}\in\bm{\Theta}}\ell(\bm{x,\delta;\theta})$.
This method is widely used due to its good properties, under regular
conditions. The MLEs are consistent, efficient and asymptotically normally
distributed, (see (Migon et al., 2014) for a detailed discussion),

$\bm{\hat{\theta}}\sim\text{N}_{k}\left(\bm{\theta},I^{-1}(\bm{\theta})\right)\quad\mbox{for}\quad
n\to\infty,$

where $I(\bm{\theta})$ is the $k\times k$ Fisher information matrix for
$\bm{\theta}$, where $I_{ij}(\bm{\theta})$ is the $(i,j)$-th element of
$I(\bm{\theta})$ given by

$I_{ij}(\bm{\theta})=\text{E}\left[-\frac{\partial^{2}}{\partial\theta_{i}\partial\theta_{j}}\ell(\bm{\theta;x,\delta})\right].$

In our specific case due to the presence of censored observations (the
censorship is random and non-informative), we cannot calculate the Fisher
information matrix $I(\bm{\theta}).$ Therefore, an alternative approach is to
use the observed information matrix $H(\bm{\theta})$ evaluated at
$\hat{\bm{\theta}}$, i.e. $H(\hat{\bm{\theta}})$. The terms of this matrix are

$H_{ij}(\hat{\bm{\theta}})=-\frac{\partial^{2}}{\partial\theta_{i}\partial\theta_{j}}\ell(\bm{\theta;x,\delta})\bigg{|}_{\bm{\theta}=\hat{\bm{\theta}}}.$

Large sample confidence intervals (approximate) at the level $100(1-\xi)\%$,
for each parameter $\theta_{i}$, $i=1,\,\ldots,\,k$, can be calculated as

$\hat{\theta}_{i}\pm Z_{\frac{\xi}{2}}\sqrt{H^{-1}_{ii}(\hat{\bm{\theta}})},$

where $Z_{\xi/2}$ denotes the $\left(\xi/2\right)$-th quantile of a standard
normal distribution. The derivation of MLEs for each of the simulation schemes
is presented in the section 5.

To perform the simulations, we consider the R language (R Core Team, 2019).
This language provides a free environment for statistical computing, which
renders many applications, such as in machine learning, text mining, and
graphs. The R language can also be compiled and run on various platforms,
being widely accepted and used in the academic community and industry. For
high-quality materials for R practitioners, see (Robert et al., 2010; Suess
and Trumbo, 2010; Field et al., 2012; James et al., 2013; Schumacker and
Tomek, 2013).

##  Sampling from random values

In this section, we revise two standard procedures to sample from a target
distribution. The procedures are simple and intuitive and can be used for most
of the probability density functions introduced in the literature.

###  The inverse transformation method

The currently most popular and well-known method used for simulating samples
from known distributions is the inverse transformation method (also known as
inverse transform sampling) (Ross et al., 2006)).

The uniform distribution in $(0,1)$ plays a central role during the
simulations procedures as, in general, the procedures are based on this
distribution. There are many methods for simulating pseudo-random samples from
uniform distributions. The congruential method and the feedback shift register
method are the most popular ones, being other procedures combinations of these
algorithms (Niederreiter, 1992). Here, we consider the standard procedure
implemented in R to generate samples from the uniform distribution.

Let $u$ be a realization of a random variable $U$ with distribution $U(0,1)$.
From the inverse transformation method, we have that $X=F^{-1}(u)$, where
$F(.)$ is a invertible distribution function, so that
$F_{X}^{-1}(u)=\inf\\{x;F(x)\geq u,u\in(0,1)\\}.$ Hence, the method involves
the computation of the quantile function of the target distribution.

The algorithm describing this methodology takes the following steps:

Algorithm 1 Inverse transformed method

1: Define $n$, and $\bm{\theta}$,

2: for $i=1$ to $n$ do

3: $u_{i}\leftarrow U(0,1)$

4: $x_{i}\leftarrow F^{-1}(u_{i},\bm{\theta})$

5: end for

To illustrate the above procedure’s applicability, we employ such a procedure
in the Weibull distribution to generate pseudo-random samples. The
calculations are performed using the language R.

###### Example 3.1

Let $X\sim\text{Weibull}(\beta,\alpha)$, to apply the inverse transform method
above to generate the pseudo-random samples we consider the cumulative density
function given by

$F(x;\beta,\alpha)=1-\exp\left\\{-\beta x^{\alpha}\right\\},\quad x>0,$

and after some algebraic manipulations, the quantile function is given by

$x_{u}=\left(\frac{\left(-\log(1-u)\right)}{\beta}\right)^{\frac{1}{\alpha}}.$
(3)

Using the R software, we can easily generate the pseudo-random sample by
considering the following algorithm:

⬇

1n=10000; beta <- 2.5; alpha <- 1.5 #Define the parameters

2rweibull<-function(n,alpha,beta) {

3 U<-runif(n,0,1)

4 t<-((-log(1-U)/beta)^(1/alpha))

5 return(t)

6}

7x<-rweibull(n, alpha, beta)

The histogram of the obtained samples is shown in Figure 1.

Figure 1: Histogram of 10,000 samples generated by the inverse transformation
method applied in the Weibull distribution. The red line shows the theoretical
probability density function shown in equation (1).

The function above will be used in all the programs discussed throughout this
paper. Hence, it should be used before running the codes.

###  Random samples for mixture models

Mixture models have been the subject of intense research due to their wide
applicability and ability to describe various phenomena whose analysis and
adjustment would be almost impossible to perform with the assumption of a
single distribution (Dellaportas et al., 1997). The mixture models were
designed to describe populations composed of sub-populations where each one
has its characteristic distribution. The mixture model’s shape is a sum of
each of the $k$ components called the mixture component, each one multiplied
by its corresponding proportion within the population. The mixing model is
expressed as follows:

$f(x)=\sum_{j=1}^{k}p_{j}f(x|\mbox{\boldmath$\theta$}_{j}),$ (4)

where the probabilities $p_{j}$, $j=1,\ldots,k$ with $\sum_{i=1}^{k}p_{i}=1,$
corresponds to the proportion of each component of the population. The
probability density $f(x)$, expressed in equation (4) with the vector
$\mbox{\boldmath$\theta$}_{j}$ of parameters, referring to the $j-$th
component of the mixture, is said to have a $k$-finite mixture density.

For the case of $k$ mixture components, we have that
$f_{1}(X),\ldots,f_{k}(X)$ are the $k$ densities of each factor with
$p_{1},\ldots,p_{k}$ proportions, respectively. The algorithm sorts these
proportions in an increasing way to generate samples from $f(x)$, so that it
is possible to identify the sample origin of the respective component,
supposing that $p_{1}<\ldots<p_{n}$. As follows, we describe the algorithm to
generate samples from a mixture model.

Algorithm 2 Simulation of random samples for mixture models

1: Define $n,k$, and $\bm{\theta}_{j}$, $p_{j}$, $j=1,\ldots k$,

2: $t$ vector

3: $p_{0}=0$

4: for $i=1$ to $n$ do

5: $u_{i}\leftarrow U(0,1)$

6: for $j=1$ to $k$ do

7: if $(p_{j-1}<u_{i}\leq p_{j})$ then

8: $t_{i}\leftarrow F_{j}^{-1}(u_{i},\bm{\theta}_{j})$

9: end if

10: end for

11: end for

12: print t

The algorithm presented above is straightforward to be applied. For instance,
to generate data for mixture models with two components $k=2$, where each
component is a Weibull, we have the following code in R.

###### Example 3.2

Let $X_{1}\sim\text{Weibull}(\beta_{1},\alpha_{1})$ and
$X_{2}\sim\text{Weibull}(\beta_{2},\alpha_{2})$ with probabilty density
function given by

$f(x|\bm{\alpha},\bm{\beta})=pf_{1}(x|\beta_{j},\alpha_{j})+(1-p)f_{2}(x|\beta,\alpha),$

where $f_{j}(x|\beta_{j},\alpha_{j})$ has density
$\text{Weibull}(\beta_{j},\alpha_{j})$ for $j=1,2$ and $p$ is the proportion
related to the first subgroup. The R codes to generate such sample given the
parameters is given as follows. In Figure 2 we show the simulatio results.

⬇

1alpha1=5;

2beta1=2.5

3alpha2=50;

4beta2=5

5n=10000;

6p1=0.8;

7p2=1-p1

8t<-c()

9x1 <- rweibull(n,alpha1,beta1)

10x2 <- rweibull(n,alpha2,beta2)

11u<-runif(n,0,1)

12for (i in 1:n) {

13 if (u[i]<=p1){

14 t[i] <- x1[i]

15 }

16 else{

17 t[i] <- x2[i]

18 }

19}

Figure 2: Histogram of 10,000 samples from a mixture of two Weibull
distributions, generated by the inverse transformation method.

###  When the inverse transformation method fails

The Markov chain Monte Carlo (MCMC) methods arise from simulating complex
distributions where the inverse transformation method fails. The method
appeared initially from applications in physic to simulate particles’ behavior
and, later, in other branches of science, including statistics. This
subsection describes the Metropolis-Hastings algorithm, which is one of the
MCMC methods, and shows how it can be used to generate random samples from a
given probability distribution.

The MH algorithm was proposed in 1953 (Metropolis et al., 1953) and further
generalized (Hastings, 1970), having the Gibbs sampler (Gelman et al., 1992)
as a particular case. This algorithm played a fundamental role in Bayesian
statistics (Chib and Greenberg, 1995). These works applied the approach both
in the context of the simulation of absolutely continuous target density and
the simulation of discrete and continuous mixed distributions. Here we provide
a brief review of the concepts needed to understand the MH algorithm. This
presentation is not exhaustive, and a detailed discussion of MCMC methods is
available in (Gamerman and Lopes, 2006) and (Robert et al., 2010) and the
references therein.

An important characteristic of the algorithm is that we can draw samples from
any probability distribution $F(\cdot)$ if we know the form of $f(x)$, even
without its normalized constant. Therefore, this algorithm is hand to generate
pseudo-random samples even when the PDF has a complex form. To obtain this
sample, we need to consider a proposal distribution of
$\operatorname{q}(\cdot)$. This proposed distribution also depends on the
value that were generated in the previous state, $x_{n-1}$, that is,
$q(\cdot|x_{n-1})$. The initial value used to initiate the algorithm does play
an important role as the Markov chain always converges to the same values. The
chain converges to the same probability distribution, usually called the
equilibrium distribution. However, note that the Markov chain must fulfill
certain properties to ensure convergence in an equilibrium distribution, such
as irreducible and reversible (Hoel et al., 1986).

The Metropolis-Hastings algorithm is summarized next.

Algorithm 3 Metropolis-Hastings algorithm

1: Start with an initial value $x^{(1)}$ and set the iteration counter $j=1;$

2: Generate a random value $x^{(*)}$ from the proposal distribution
$Q(\cdot)$;

3: Evaluate the acceptance probability
$\alpha\left(x^{(j)},x^{(*)}\right)=\min\left(1,\frac{f\left(x^{(*)}|\bm{x}\right)}{f\left(x^{(j)}|\bm{x}\right)}\frac{\operatorname{q}\left(x^{(j)},x^{(*)},b\right)}{\operatorname{q}\left(x^{(*)},x^{(j)},b\right)}\right),$
where $f(\cdot)$ is the probability distribution of interest. Generate a
random value $u$ from an independent uniform in $(0,1)$;

4: If $\alpha\left(x^{(j)},x^{{(*)}}\right)\geq u(0,1)$ then
$x^{(j+1)}=x^{(*)}$, otherwise $x^{(j+1)}=x^{(j)}$;

5: Change the counter from $j$ to $j+1$ and return to step 2 until convergence
is reached.

The proposal distribution should be selected carefully. For instance, if we
want to sample from the Weibull distribution where $x>0$, using the Beta
distribution with $0\leq x\leq 1$ is not a good choice as they have different
domains. In this case, a common approach is to consider a truncated normal
distribution with a range similar to the target distribution domain. Further,
the pseudo-random sample generated from the algorithm above can be considered
as a sample of $f(\cdot)$.

###### Example 3.3

The Weibull distribution can be written as
$f(x;\beta,\alpha)=C(\beta,\alpha)x^{\alpha-1}\exp\left\\{-\beta
x^{\alpha}\right\\}$ where $C$ is a normalizing function. Suppose that we drop
the power parameter $\alpha$ in $e^{-\beta x^{\alpha}}$, remove $x^{-1}$ and
include the possibility of the presence of a lower bound for $x$, refereed as
$x_{\min}$. Thus, we obtain

$f(x\,|\alpha,\beta,x_{\min})=C(\alpha,\beta,x_{\min})x^{-\alpha}e^{-\beta
x},$ (5)

where $\alpha>0$ and $\beta>0$ are the unknown shape and scale parameters,
$C(\cdot)$ is the normalized constant given by

$C(\alpha,\beta,x_{\min})=\frac{\beta^{1-\alpha}}{\Gamma(1-\alpha,\beta
x_{\min})},$ (6)

where

$\Gamma(s,x)=\int_{x}^{\infty}t^{s-1}e^{-t}dt,$ (7)

is the upper incomplete gamma function. This distributions is known as the
power-law distribution with cutoff (Clauset et al., 2009) and has an more
complicate structure than its related power-law distribution. Note that in the
limit when $\beta\rightarrow 0$ the cutoff term is

$\lim_{\beta\to 0}=\frac{\beta^{1-\alpha}}{\Gamma(1-\alpha,\beta
x_{\min})}=\frac{\alpha-1}{x_{\min}}\left(\frac{1}{x_{\min}}\right)^{-\alpha}.$

From the expression above and as the exponential term tends to unity when
$\beta\rightarrow 0$, we obtain the standard power-law distribution given by

$f(x\,|\alpha,x_{\min})=\frac{\alpha-1}{x_{\min}}\left(\frac{x}{x_{\min}}\right)^{-\alpha}.$
(8)

The cumulative distribution function (CDF) of $f(x\,|\alpha,\beta,x_{\min})$
has the form,

$\displaystyle F(x\,|\alpha,\beta,x_{\min})$
$\displaystyle=\frac{\beta^{1-\alpha}}{\Gamma(1-\alpha,\beta
x_{\min})}\int_{x_{\min}}^{x}y^{-\alpha}e^{-\beta y}dy.$

Since the CDF does not have a closed-form expression, we cannot obtain the
expression of the quantile function of the power-law model in a closed-form.
While in most situations, we can obtain pseudo-random samples using a root
finder with $u_{i}-F(x_{i})=0$, the obtained samples for the power-law
distribution with cutoff very different from the expected ones. The problem
may occur due to the integrals involved in the CDF. To overcome this problem,
we consider the Metropolis-Hasting algorithm. A function in R to sample from a
power-law distribution can be seen as follows.

⬇

1library(truncnorm)

2library(TruncatedNormal)

3

4n<-10000; xm<-1; alpha<-1.5 ; beta<-0.5 ;

5rplc<-function(n,xm,alpha,beta,burnin=1000,thin=50,se=0.5){

6 R<-n*thin+burnin #Number of replication

7 #Density without the normalized constant

8 den <- function (x) {-alpha*log(x)-beta*x}

9 xv<-length(R+1);

10 xv[1]<-xm+1;

11 for(i in 1:R) {

12 prop1<- rtnorm(1,xv[i],se,lower=xm, upper=Inf)

13 dens1<-dtnorm(xv[i], mean=prop1, sd=se, lower=xm, upper=Inf,

14 log = TRUE)

15 dens2<-dtnorm(prop1, mean=xv[i], sd=se,lower=xm,

16 upper=Inf, log = TRUE)

17 ratio1<-den(prop1)-den(xv[i]) +dens1-dens2

18 has<-min(1,exp(ratio1)); u1<-runif(1)

19 if (u1<has & is.double(has)) {

20 xv[i+1]<-prop1

21 }

22 else {

23 xv[i+1]<- xv[i]

24 }

25 }

26 x<-xv[seq(burnin,R,thin)]

27 return(x)

28}

29x<-rplc(n,xm,alpha,beta)

The function rplc can be used to obtain a pseudo-random sample from the power-
law with a cutoff (PLC) distribution. In this case, some additional parameters
need to be set: (i) the burning, (ii) the thin which is used to decrease the
auto-correlation between the elements of the sample, and (iii) the se, i.e.,
the standard error of the truncated normal distribution chosen as $q(\cdot)$.
This distribution was assumed due to its flexibility to set a lower bound
$x_{min}$, hence provides values in the same domain of the PLC distribution.
For the additional parameters, we set default values for the burning $=1000$,
thin $=50$, and se $=0.5$. Hence, if different values are not presented, the
function will consider these values to generate the pseudo-random samples.
Figure 3 shows the histogram of 10,000 samples generated by this function and
the respective theoretical curve, given by equation (5).

Figure 3: Histogram of 10,000 samples from the power-law distribution in
equation (5), generated by the Metropolis-Hastings algorithm.

To verify if the estimates of the PLC parameters are closer to the true
values, we derived the MLEs for the parameter of the PLC distribution, which
can be seen in Appendix A.

##  Sampling right-censored data

Sampling pseudo-random data with censoring is usually a challenging task. In
this section, we discuss simple procedures to generate such pseudo-random
samples in the presence of type I, II, and random censoring. These censoring
schemes are the most common in practice. In the end, we also present some
references to generate other types of censoring.

###  Type II censorship

To generate the type II censorship, we define the total number of observations
$n$ and the number $m$ of observations censored. First, we generate a sample
of size $n$ according to a given distribution. We sort these observations in
increasing order and associate the first $n-m$ values to the observed data.
The remainder $m$ observations receives the value in position $m+1$. Since
these observations did not experience the event of interest, we assume that
their survival time is equal to the largest one observed. The algorithm is
described as follows.

Algorithm 4 Simulation With Type II censorship

1: Define $n$, $m$ and $\bm{\theta}$,

2: for $i=1$ to $n-m$ do

3: $u_{i}\leftarrow U(0,1)$

4: $x_{i}\leftarrow F^{-1}(u_{i},\bm{\theta})$

5: $\delta_{i}\leftarrow 1$

6: end for

7: $t\leftarrow\text{sort}(x)$

8: for $i=n-m+1$ to $n$ do

9: $t_{i}\leftarrow t_{n-m}$

10: $\delta_{i}\leftarrow 0$

11: end for

12: print $t$, $\delta$

The algorithm above is general and can be used for any probability
distribution from where one can obtain a pseudo-random sample. Note that we
use the inverse transformation method to obtain such a sample. However, the
same procedure can be used, assuming other sampling schemes. To illustrate the
proposed approach, let us consider an example with the Weibull distribution.

###### Example 4.1

Let $X\sim\text{Weibull}(\beta,\alpha)$. We can generate samples with type II
censoring by fixing the number of elements to be censored. We will generate
$n=50$ observations from the Weibull distribution and set the number of
censoring $m=15$. The code in R is given as follows.

⬇

1######Weibull R data simulation With Type II censorship#####

2n=50 #Sample size

3m=15 #Number of censored data

4beta <- 2.5 ##Parameter value

5alpha <- 1.5 ##Parameter value

6x <- rweibull(n,alpha,beta)

7t <- sort(x)

8r<-n-m

9t[(r+1):n]=t[r]

10delta<-ifelse(seq(1:n)<(n-m+1),1,0) #Vector indication censoring

An example of the output obtained from the code above in R is given by:

    > print(t)
     [1] 0.01 0.01 0.03 0.09 0.10 0.10 0.15 0.15 0.17 0.17 0.20 0.22 0.22 0.25 0.31
    [16] 0.32 0.34 0.34 0.35 0.36 0.39 0.40 0.41 0.44 0.46 0.46 0.47 0.49 0.53 0.53
    [31] 0.54 0.56 0.56 0.62 0.67 0.67 0.67 0.67 0.67 0.67 0.67 0.67 0.67 0.67 0.67
    [46] 0.67 0.67 0.67 0.67 0.67
    > print(delta)
     [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0
    [39] 0 0 0 0 0 0 0 0 0 0 0 0

We can see that 15 observations are censored and receive the largest lifetime
among the observed lifetimes. Also, the vector _delta_ stores the
classification of each sample, where the value 0 is attributed to censored
data.

###  Type I censorship

The type I censoring data is commonly observed in different applications
related to medical and engineering studies. Here, we have a censoring time
that the experiment will be ended, while in type II censorship the study ends
only after ”$n-m$” failures. For instance, suppose that the reliability of a
component will be tested in a controlled study. Since many components are
subject to high reliability, the experiment could take many years to be finish
and the product would be outdated. To overcome this problem, we define the
final time of the experiment, thus, given the $n$ observations generated from
the distribution under study, we define a time $t_{c}$ so that every component
with survival time greater than the fixed time will correspond to a censored
observation. To generate this type of censorship, consider the following
algorithm.

Algorithm 5 Simulation With Type I censorship

1: Define $n$, $t_{c}$ and $\bm{\theta}$,

2: for $i=1$ to $n$ do

3: $u_{i}\leftarrow U(0,1)$

4: $x_{i}\leftarrow F^{-1}(u_{i},\bm{\theta})$

5: if $x_{i}<t_{c}$ then

6: $t_{i}\leftarrow x_{i}$

7: $\delta_{i}\leftarrow 1$

8: else

9: $t_{i}\leftarrow t_{c}$

10: $\delta_{i}\leftarrow 0$

11: end if

12: end for

13: print $t$, $\delta$

As it can be noted, we defined $\delta$ as the indicator of censoring, i.e.,
if $\delta=1$ the component has observed the event of interest, whereas
$\delta=0$ indicates that the component is censored.

###### Example 4.2

Let us consider again that $X\sim\text{Weibull}(\beta,\alpha)$, to simulate
data with type I censorship. We set the parameter values and the time that the
experiment $t_{c}$ will end. Assuming that $\alpha=1.5$, $\beta=2.5$ and
$t_{c}=1.5$ the codes to generate the censored data are given as follows.

⬇

1#####Weibull R data simulation With Type I

2n=50 #Sample size

3tc<- 0.51 ##Censored time

4beta <- 2.5 ##Parameter value

5alpha <- 1.5 ##Parameter value

6t<-rweibull(n, alpha, beta)

7delta<-rep(1,n)

8for (i in 1:n)

9 if(t[i]>=tc){

10 t[i]<-tc

11 delta[i]<-0 }

The obtained output in R is given by

    > print(t)
     [1] 0.51 0.34 0.51 0.41 0.15 0.09 0.15 0.34 0.01 0.51 0.51 0.51 0.51 0.36 0.35
    [16] 0.46 0.51 0.51 0.46 0.25 0.20 0.10 0.51 0.51 0.51 0.17 0.31 0.51 0.39 0.47
    [31] 0.51 0.10 0.51 0.51 0.01 0.32 0.03 0.51 0.44 0.22 0.51 0.22 0.51 0.51 0.51
    [46] 0.49 0.40 0.51 0.51 0.17
    > print(delta)
     [1] 0 1 0 1 1 1 1 1 1 0 0 0 0 1 1 1 0 0 1 1 1 1 0 0 0 1 1 0 1 1 0 1 0 0 1 1 1 0
    [39] 1 1 0 1 0 0 0 1 1 0 0 1

An important concern is how do we set $t_{c}$ to obtain a desirable level of
censored data. For instance, what is the ideal $t_{c}$ to obtain $20\%$ of
censored observations? Although in applications the $t_{c}$ is usually
defined, here, we have to find the value of $t_{c}$ to generate the pseudo-
random samples. To achieve that, we can consider the inverse transformed
method since,

$F(t_{c};\bm{\theta})=1-\pi\quad\Rightarrow\quad
t_{c}=F^{-1}(\pi;\bm{\theta}).$

Hence, if we want $40\%$ censorship, it is necessary to find the $t_{c}$ in
the quantile function at the point $0.6$, i.e.,

$\int_{0}^{t_{c}}f_{X}(x)dx=1-\pi,$

where $\pi$ is the percentage of censorship. Then if $x>t_{c}$ is censorship
$x\leq t_{c}$ is complete then we expect $0.6$ to be complete and $0.4$
censored. For example, in the case where your data is generated from the
Weibull distribution, we will have the following,

$\int_{0}^{t_{c}}\alpha\beta x^{\alpha-1}\exp\\{-\beta x^{\alpha}\\}dx=1-\pi.$
(9)

The quantile function for our scheme takes the form

$t_{c}=\left(\frac{\left(-\log(\pi)\right)}{\beta}\right)^{\frac{1}{\alpha}}.$

Hence, for a value of $\pi=0.4,\beta=2.5,\alpha=1.5$, we would have that
$t_{c}=0.5121$. In fact, running the code of Example 4.2, with n=50000 we
obtain $\hat{\pi}=0.3969$.

###  Random censoring

In the simulation method with random censorship for survival data, we need to
generate two distributions. The first one controls the survival times without
censorship, whereas the second distribution controls the censorship mechanism.
Based on these elements, the censored event is the one in which the uncensored
survival time generated in the first distribution is greater than the censored
time generated in the second distribution.

Specifically, consider a non-negative random variable $X$ representing our
first distribution that generates the times until the event occurs. The
survival function $\mathcal{S}_{X}(x)$ is defined from the cumulative
distribution function $F_{X}(x)$ and offers the probability of survival after
time $x$. That is, $F_{X}(x)=P(X\leq x)$ and,
$\mathcal{S}_{X}(x)=1-P(X\leq~{}x)=P(X>x)$ where $f_{X}(x)$ is the probability
density function of $X$. For the second distribution that controls the
censorship mechanism $C$, the choice is at the discretion of the researcher.
However, it is recommended to consider a continuous distribution. Here, we
assume that the random variable follows the $\text{Weibull}(\beta,\alpha)$ and
$C$ is uniformly distributed, $C\sim U(0,\lambda)$. We also assume that these
two distributions are independent $(X\perp U)$. Therefore, our censored
samples have the following scheme. For each $x$-th observation, we have
$(T_{x},\delta_{x})$ with $T_{x}=\min(X_{x},C_{x})$ and
$\delta_{x}=\mathrm{I}(X_{x}\leq C_{x})$. Where $\delta=1$ when the lifetime
is observed or $\delta=0$, when it is censored. The following algorithm
describes the random censoring.

Algorithm 6 Simulation With Random censoring

1: Define $n$, $\lambda$ and $\bm{\theta}$,

2: for $i=1$ to $n$ do

3: $u_{i}\leftarrow U(0,1)$

4: $x_{i}\leftarrow F^{-1}(u_{i},\bm{\theta})$

5: $c_{i}\leftarrow U(0,\lambda)$

6: if $x_{i}<c_{i}$ then

7: $t_{i}\leftarrow x_{i}$

8: $\delta_{i}\leftarrow 1$

9: else

10: $t_{i}\leftarrow c_{i}$

11: $\delta_{i}\leftarrow 0$

12: end if

13: end for

14: print $t$, $\delta$

###### Example 4.3

Assuming that the distribution that controls censorship $C\sim U(0,\lambda)$,
we generate data with random censorship for the time generated from the
$X\sim\text{Weibull}(\beta,\alpha)$ distribution. The code in R is described
as follows.

⬇

1n=50; beta <- 2.5; alpha <- 1.5 #Define the parameters

2lambda<-1.22

3t<-c()

4y<-rweibull(n, alpha, beta)

5delta<-rep(1,n)

6c<-runif(n,0,lambda)

7for (i in 1:n) {

8 if (y[i]<=c[i]) {

9 t[i]<-y[i];

10 delta[i]<-1

11 }

12 else {

13 t[i]<-c[i];

14 delta[i]<-0

15 }

16 }

An example of output is given as follows.

    > print(t)
     [1] 0.56 0.34 0.53 0.41 0.15 0.00 0.15 0.27 0.01 0.53 0.25 0.67 0.28 0.36 0.35
    [16] 0.46 1.19 0.56 0.46 0.20 0.20 0.10 0.47 0.57 0.65 0.17 0.31 0.07 0.39 0.47
    [31] 1.00 0.10 0.37 0.48 0.01 0.21 0.03 0.44 0.15 0.22 0.73 0.22 0.43 0.91 0.54
    [46] 0.22 0.40 0.19 0.74 0.17
    > print(delta)
     [1] 1 1 1 1 1 0 1 0 1 1 0 1 0 1 1 1 1 1 1 0 1 1 0 0 0 1 1 0 1 1 0 1 0 0 1 0 1 0
    [39] 0 1 0 1 0 1 1 0 1 0 1 1

To find the desired censorship percentage in the random censorship scheme, the
following methodology is proposed in Martinez et al. (2016). It is assumed
independence between the random variables $C\sim U(0,\lambda)$ and
$X\sim\text{Weibull}(\beta,\alpha)$ that represents for each observation
$x$-th the associated censorship and the times until the event occurred
respectively, with $\theta$ being the desired percentage of censorship. Note
that in the case of censored data, we have the following situation,
$(C_{x},0)\Leftrightarrow T_{x}=C_{x}\Leftrightarrow C_{x}<X_{x}$. Hence, by
definition, $P(C_{x}<X_{x})$, for $C_{x}\perp X_{x}$, is given by

$\displaystyle P(C<X)$
$\displaystyle=\int_{0}^{\infty}\int_{0}^{x}f_{X,C}(x,c)dcdx=\int_{0}^{\infty}\int_{0}^{x}f_{X}(x)f_{C}(c)dcdx\int_{0}^{\infty}F_{C}(x)f_{X}(x)dx=\theta,$
(10)

In the present paper, we consider the Weibull distribution for $X$ and Uniform
distribution for $C$. Therefore,

$\displaystyle P(C<X)=\int_{0}^{\infty}\tfrac{x}{\lambda}\alpha\beta
x^{\alpha-1}\exp\left\\{-\beta
x^{\alpha}\right\\}dx=\int_{0}^{\infty}\tfrac{\alpha\beta}{\lambda}x^{\alpha}\exp\left\\{-\beta
x^{\alpha}\right\\}dx=\pi.$

The solution of the equation above can be easily obtained after some algebraic
manipulation that leads to

$\lambda=\frac{\beta^{-\frac{1}{\alpha}}}{\alpha\pi}\Gamma\left(\frac{1}{\alpha}\right).$
(11)

The solution of equation (11) provides us the necessary value to achieve a
proportion of censoring $\pi$. For instance, recalling Example 4.3, assuming
that the same parameter values $\alpha=1.5$ and $\beta=2.5$ but seating
$\pi=0.4$, we have that $\lambda=1.22$, which is the same value used in the
cited example.

The code in R for this example is given as follows.

⬇

1n=100000; beta <- 2.5; alpha <- 1.5 #Define the parameters

2theta<-0.4

3lambda<-(beta^(-1/alpha))*gamma(1/alpha)/(alpha*theta)

4t<-c()

5y<-rweibull(n, alpha, beta)

6delta<-rep(1,n)

7c<-runif(n,0,lambda)

8for (i in 1:n) {

9 if (y[i]<=c[i]) {

10 t[i]<-y[i];

11 delta[i]<-1

12 }

13 else {

14 t[i]<-c[i];

15 delta[i]<-0

16 }

17 }

18print(1-(sum(delta)/n))

The output is given by,

    > 1-(sum(delta)/n) ##Estimated theta
    [1] 0.3894

Note that we have considered a large sample size since, as the sample
increases, the censoring percentage tends to be true. If we sample many
samples and take the average, it is expected to achieve a mean of $\pi$. A
major limitation of this approach is that in many models, the PDF $f_{X}(x)$
has a complex structure, and the solution cannot be obtained easily. Although
we can use numerical methods to find such a solution, we have faced failure to
convergence the method or values that do not return the desirable $\pi$
proportion of censored data. To overcome this problem, we have considered a
simple grid procedure to obtain such value.

###  Sampling with cure fraction

In medical and public health studies, it is common to observe a proportion of
individuals not susceptible to the event of interest. These individuals are
usually called immune or cured, and the fraction of such individuals is called
cure fraction. The survival models that accommodate this characteristic are
referred to as cure rate models or long-term models (see, for instance, Maller
and Zhou (1996)).

Rodrigues et al. (Rodrigues et al., 2009) proposed a unified version of long-
term survival models. In this case, let $M$ be the number of competing causes
related to the occurrence of an event of interest. Suppose that $M$ follows a
discrete distribution with a probability mass function (pmf) given by
$p_{m}=P(M=m)$, for $m=0,1,2,\ldots$. Given $M=m$, consider
$Y_{1},Y_{2},\ldots,Y_{m}$ independent identically distributed (iid)
continuous non-negative random variables that do not depend on $M$ and with a
common survival function $S(y)=P(Y\geq y)$. The random variables $Y_{k}$,
$k=1,2,\ldots,m$, denote the time-to-failure of a individual due to the $k$-th
competing cause. Thus, the observable lifetime is defined as
$X=\min\\{Y_{1},Y_{2},\ldots,Y_{M}\\}$ for $M\geq 1$, with
$P(Y_{0}=+\infty)=1$, which leads to a proportion $p_{0}$ of the population
not susceptible to the failure. Hence, the long-term survival function of the
random variable $X$ is given by

$\displaystyle S_{\rm pop}(x)$ $\displaystyle=$ $\displaystyle P\left(X\geq
x\right)$ (12) $\displaystyle=$ $\displaystyle
P\left(M=0\right)+P\left(Y_{1}>x,Y_{2}>x,\ldots,Y_{M}>x,M\geq 1\right)$
$\displaystyle=$ $\displaystyle
p_{0}+\sum_{m=1}^{\infty}p_{m}P\left(Y_{1}>x,Y_{2}>x,\ldots,Y_{M}>x\mid
M=m\right)$ $\displaystyle=$
$\displaystyle\sum_{m=0}^{\infty}p_{m}\left[S(x)\right]^{m}$ $\displaystyle=$
$\displaystyle G_{M}\left(S(x)\right),$

where $S(x)$ is a proper survival function associated with the individuals at
risk in the population and $G_{M}\left(\cdot\right)$ is the probability
generating function (PGF) of the random variable $M$, whose sum is guaranteed
to converge because $S(t)\in[0,1]$ (see Feller (2008)).

Suppose that the number of competing causes follows a Binomial negative
distribution in equation (12), with parameters $\kappa$ and $\gamma$, for
$\gamma>0$ and $\kappa\gamma\geq-1$, following the parameterization of
Piegorsch (1990). In this case, the long-term survival model is given by

$\displaystyle
S_{pop}(x)=\left(1+\kappa\gamma(1-S(x))\right)^{-1/\kappa},\quad x>0.$ (13)

The Negative Binomial distribution includes the Bernoulli ($\kappa=-1$ and
$m=1$) and Poisson ($\kappa=0$) distributions as particular cases. Hence,
equation (13) reduces to the mixture cure rate model, proposed by Berkson and
Gage (1952), and to the promotion time cure model, introduced by Tsodikov et
al. (1996), respectively. Here, we focus on the mixture cure rate model, since
most of the cure rate models can be written as

$\displaystyle S_{pop}(x)=p+(1-p)S(x),\quad x>0\quad\mbox{and}\quad 0\leq
p\leq 1.$ (14)

An algorithm to generate a random sample for the model above is shown as
follows.

Algorithm 7 Sampling With Cure fraction

1: Define $n$, $\beta$, $\alpha$, $p$ and $\lambda$,

2: for $i=1$ to $n$ do

3: $b_{i}\leftarrow\text{Ber}(n,p)$

4: if $b_{i}=1$ then

5: $y_{i}\leftarrow\infty$

6: else

7: $u_{i}\leftarrow U(0,1)$

8: $y_{i}\leftarrow F^{-1}(u_{i},\bm{\theta})$

9: end if

10: end for

11: for $i=1$ to $n$ do

12: $c_{i}\leftarrow U(0,\lambda)$

13: if $x_{i}<c_{i}$ then

14: $t_{i}\leftarrow x_{i}$

15: $\delta_{i}\leftarrow 1$

16: else

17: $t_{i}\leftarrow c_{i}$

18: $\delta_{i}\leftarrow 0$

19: end if

20: end for

21: print $t$, $\delta$

###### Example 4.4

Let $X\sim\text{Weibull}(\beta,\alpha)$ and suppose that the population is
composed of two groups: (i) cured, with a cure fraction $p$, and (ii) non-
cured, with a cure fraction $1-p$ and survival function
$S_{0}(x;\alpha,\beta)=\exp\\{-\beta x^{\alpha}\\}$. Thus, the survival
function of the entire population, $S_{p}(x;\alpha,\beta)$, is given by:

$\displaystyle S_{p}(x;\alpha,\beta)$ $\displaystyle=p+(1-p)\exp\\{-\beta
x^{\alpha}\\}.$ (15)

The pseudo-random sample with cure fraction using R can be generated using the
following code.

⬇

1n=50; beta <- 2.5; alpha <- 1.5; p<-0.3 #Define the parameters

2lambda<-1.22

3

4require(Rlab)

5##Generation of pseudo-random samples with long-term survival

6b<-rbern(n,p)

7t<-c(); y<-c()

8for (i in 1:n) {

9 if (b[i]==1) {

10 y[i]<-1.000e+54

11 }

12 else{

13 y[i]<-rweibull(1, alpha, beta)

14 }

15}

16##Generation of pseudo-random samples with random censoring

17delta<-c()

18cax<-runif(n,0,lambda)

19for (i in 1:n) {

20 if (y[i]<=cax[i]) {t[i]<-y[i] ;delta[i]<-1}

21 if (y[i]>cax[i]) {t[i]<-cax[i] ;delta[i]<-0}

22}

The obtained output in R is given by

    > print(t)
     [1] 0.48 0.26 0.21 0.29 0.44 0.01 0.27 0.22 0.46 0.41 0.97 0.62 0.22 0.21 0.19
    [16] 0.22 0.52 0.37 0.32 0.53 1.14 0.31 0.20 0.83 0.96 0.56 0.18 0.67 0.04 0.52
    [31] 1.13 0.11 0.49 0.35 0.40 0.45 0.27 1.19 0.26 0.08 1.02 0.07 1.09 1.22 0.63
    [46] 0.77 0.12 0.79 0.63 0.28
    > print(delta)
     [1] 0 1 0 0 0 1 0 1 1 1 0 0 0 1 0 1 0 1 0 1 0 0 0 0 0 1 1 0 0 0 0 0 0 0 1 1 1 0
    [39] 1 1 0 0 0 0 1 1 1 0 0 1 1

As can be seen in the results above, the proportion of censoring was $0.6$,
which is different from $0.4$, that would be achieved by selecting
$\lambda=1.22$. Since we are dealing with a new model with an additional
parameter, the results need to be computed using the new PDF. On the other
hand, in most cases, we will not be able to find closed-form solutions as
obtained in (11) and the results from numerical integration may not return the
desirable proportion of censoring. To overcome this problem we can consider a
grid search approach to find the value of $\lambda$, as illustrated in the
algorithm below.

Algorithm 8 Fiding the desariable censoring

1: Define $n$, $\bm{\theta}$, $\theta_{I}$, $\lambda$,

2: repeat

3: $\lambda=\lambda+\epsilon$

4: Call Algorithm 6

5: $\hat{\theta}=1-\frac{\sum_{i=1}^{n}(\delta_{i})}{n}$

6: until $\theta_{I}\leq\hat{\theta}$

7: print $\lambda$

In the previous algorithm, it is necessary to set an initial value for
$\lambda$, the ideal proportion of censoring $\theta_{I}$, and call Algorithm
6 during the iterative procedure. This search is straightforward, and the
computational time required to find the desired value of $\theta$ may increase
depending on the initial value of $\lambda$ and the $\epsilon$ that is used to
increase the value of $\lambda$. We can set $\lambda$ as the minimum value
from a simulated sample, which provides a better initial value than 0. The
epsilon value may depend on the range of the data since, in most applications,
we are dealing with time we usually consider $\epsilon=0.01$. Finally, to
obtain a good precision, $n$ must be very large — in practice, we usually set
$n=10,000$. It is essential to point out that although we have presented this
approach for long-term survival models, the same approach can be used for
random censoring and type I censored data.

⬇

1library(Rlab)

2n=10000; beta <- 2.5; alpha <- 1.5; p<-0.3 #Define the parameters

3lambda<-min(rweibull(n, alpha, beta))

4epsilon<-0.01

5thetai<-0.6

6theta<-1 #Initial value of theta

7

8##Creting pseudo-random samples with long-term survival

9while(thetai<=theta) {

10 lambda<-lambda+epsilon

11 b<-rbern(n,p)

12 t<-c(); y<-c()

13 for (i in 1:n) {

14 if (b[i]==1) {y[i]<-1.000e+54}

15 if (b[i]==0) {y[i]<-rweibull(1, alpha, beta)}

16 }

17 ##Creting pseudo-random samples with random censoring

18 delta<-c()

19 cax<-runif(n,0,lambda)

20 for (i in 1:n) {

21 if (y[i]<=cax[i]) {t[i]<-y[i] ;delta[i]<-1}

22 if (y[i]>cax[i]) {t[i]<-cax[i] ;delta[i]<-0}

23 }

24 theta<-1-(sum(delta)/n)

25}

The obtained output in R is given by

    > print(theta)
     [1] 0.5938
    > print(lambda)
     [1] 1.110444

##  Estimating the parameters under censoring data

Here, we revisit the MLEs for the Weibull distribution under different types
of censoring and long-term survival. The results presented are not new in the
literature and are standard in most statistical books related to survival
analysis. The aim here is to estimate the distribution parameters using the
pseudo-random samples presented in Section 4.

###  Censorship type II

Let $X_{1},\cdots,X_{n}$ be a random sample from Weibull$(\beta,\alpha)$. For
the type II censoring, the likelihood function is given by:

$L(\beta,\alpha;\bm{x})=\frac{n!}{(n-r)!}(\alpha\beta)^{r}\exp\left\\{-\beta\sum_{i=1}^{r}x_{(i)}^{\alpha}\right\\}\left(\prod_{i=1}^{r}x_{(i)}\right)^{\alpha-1}\exp\left\\{-\beta(n-r)x_{(r)}^{\alpha}\right\\}.$

The log-likelihood function without the constant term is given by:

$\ell(\beta,\alpha;\bm{x})\propto
r\left[\log(\alpha)+\log(\beta)\right]-\beta\sum_{i=1}^{r}x_{(i)}^{\alpha}+(\alpha-1)\sum_{i=1}^{r}\log(x_{(i)})-\beta(n-r)x_{(r)}^{\alpha}.$

Setting $\dfrac{d}{d\beta}\ell(\beta,\alpha;\bm{x})=0$,
$\dfrac{d}{d\alpha}\ell(\beta,\alpha;\bm{x})=0$ and after some algebraic
manipulations, we obtain the MLE $\hat{\beta}$ of $\beta$,

$\hat{\beta}=\dfrac{r}{\sum_{i=1}^{r}x_{(i)}^{\hat{\alpha}}+(n-r)x_{(r)}^{\hat{\alpha}}},$

where the MLE $\hat{\alpha}$ can be obtained as solution of equation:

$\frac{r}{\alpha}+\sum_{i=1}^{r}\log(x_{(i)})=\left(\dfrac{r}{\sum_{i=1}^{r}x_{(i)}^{\alpha}+(n-r)x_{(r)}^{\alpha}}\right)\left[\sum_{i=1}^{r}x_{(i)}^{\alpha}\log(x_{(i)})+(n-r)x_{(r)}^{\alpha}\log(x_{(r)})\right].$

The non-linear equation above has to be solved by numerical methods, such as
Newton-Raphson, to find $\hat{\alpha}$. In R, we can use the next code to
solve this equation.

⬇

1

2likealpha<-function(alpha){

3r/alpha+sum(log(w))-(r/(sum(w^alpha)+(n-r)*w[r]^alpha))*(sum(w^alpha*log(w))+

4(n-r)*(w[r]^alpha)*log(w[r]))

5}

6

7###Simplifying the vector###

8w<-t[1:r]

9

10###Computing the MLE###

11alpha<-uniroot(f=likealpha, c(0,100))$root

12beta<-r/(sum(w^alpha)+(n-r)*w[r]^alpha)

To facilitate the computation, we have kept only the complete sample in the
$w$ vector. The information related to the partial information is included
directly by w[r]. The non-linear equation was included in a function that,
combined with the uniroot function, can find the value of $\hat{\alpha}$.
Using the same data generated in Example 4.1 and running the code above, the R
output is given by

    > print(c(alpha,beta))
    [1] 1.25 1.92

Hence, the obtained estimates are $\hat{\alpha}=1.25$ and $\hat{\beta}=1.92$.
These estimates differ when compared with the true value $\alpha=1.5$ and
$\beta=2.5$. However, as the sample size increases, the values are closer to
the true value. For instance, if we consider a larger sample of $n=50,000$ and
$m=15000$, the estimates for the parameter will be $\hat{\alpha}=1.49$ and
$\hat{\beta}=2.50$ (the codes are available in the Supplemental Material).

###  Censorship type I

For type I censoring, we consider a random sample $X_{1},\cdots,X_{n}$ from a
Weibull distribution. In this case, the likelihood function is given by:

$L(\beta,\alpha;\bm{x})=\exp\left\\{-\beta(n-r)x_{c}^{\alpha}\right\\}\prod_{i=1}^{n}\left(\alpha\beta
x_{i}^{\alpha-1}\exp\left\\{-\beta
x_{i}^{\alpha}\right\\}\right)^{\delta_{i}}.$

Then, the log-likelihood function can be expressed by:

$\ell(\beta,\alpha;\bm{x})=-\beta(n-r)x_{c}^{\alpha}+r\left[\log(\alpha)+\log(\beta)\right]+(\alpha-1)\sum_{i=1}^{n}\delta_{i}\log(x_{i})-\beta\sum_{i=1}^{n}\delta_{i}x_{i}^{\alpha}.$

Considering $\dfrac{d}{d\beta}\ell(\beta,\alpha;\bm{x})=0$ and
$\dfrac{d}{d\alpha}\ell(\beta,\alpha;\bm{x})=0$, we find the likelihood
equations given by:

$\frac{r}{\beta}=(n-r)x_{c}^{\alpha}+\sum_{i=1}^{n}\delta_{i}x_{i}^{\alpha}$

and

$\frac{r}{\alpha}+\sum_{i=1}^{n}\delta_{i}\log(x_{i})=\beta(n-r)x_{c}^{\alpha}\log(x_{c})+\beta\sum_{i=1}^{n}\delta_{i}x_{i}^{\alpha}\log(x_{i}).$

Hence, we obtain the MLE of $\beta$, $\hat{\beta}$ say, given by:

$\hat{\beta}=\dfrac{r}{(n-r)x_{c}^{\hat{\alpha}}+\sum_{i=1}^{n}\delta_{i}x_{i}^{\hat{\alpha}}},$

where $\hat{\alpha}$ is the MLE of $\alpha$ obtained by solving the nonlinear
equation:

$\frac{r}{\alpha}+\sum_{i=1}^{n}\delta_{i}\log(x_{i})=\left(\dfrac{r}{(n-r)x_{c}^{\alpha}+\sum_{i=1}^{n}\delta_{i}x_{i}^{\alpha}}\right)\left[(n-r)x_{c}^{\alpha}\log(x_{c})+\sum_{i=1}^{n}\delta_{i}x_{i}^{\alpha}\log
x_{i}\right].$

The MLE for type I censoring has a similar structure as the one presented
under type II censored data. Although here $m=n-r$ is random, the MLEs can be
computed using similar techniques, as discussed in Section 5.1. The codes to
solve the equation are presented as follows.

⬇

1likealpha<-function(alpha){

2r/alpha+sum(delta*log(t))-(r/(sum(delta*(t^alpha))+

3(n-r)*tc^alpha))*((n-r)*(tc^alpha)*log(tc)+sum(delta*(t^alpha)*log(t)))

4}

5

6r<-sum(delta)

7###Computing the MLE###

8alpha<-try(uniroot(f=likealpha, c(0,100))$root)

9beta<-r/(sum(delta*t^alpha)+(n-r)*tc^alpha)

Using the same data that were generated in Example 4.2 and running the code
above, the output is given by

    > print(c(alpha,beta))
    [1] 1.26 1.96

The estimates for the parameters are $\hat{\alpha}=1.26$ and
$\hat{\beta}=1.96$. As in the previous example, as the sample size increases,
the values get closer to $\alpha=1.5$ and $\beta=2.5$. For example,
considering $n=50,000$ the estimates for the parameter would be
$\hat{\alpha}=1.50$ and $\hat{\beta}=2.52$ with $39.71\%$ of censored data.

###  Random censorship

Assuming that $X_{1},\cdots,X_{n}$ is a random sample from the Weibull
distribution, the likelihood function considering data with random censoring
can be written as:

$\displaystyle L(\beta,\alpha;\bm{x})=\prod_{i=1}^{n}\left(\alpha\beta
x_{i}^{\alpha-1}\right)^{\delta_{i}}\exp\left\\{-\beta
x_{i}^{\alpha}\right\\}.$

The log-likelihood function is given by:

$\ell(\beta,\alpha;\bm{x})=r\left[\log(\alpha)+\log(\beta)\right]+(\alpha-1)\sum_{i=1}^{n}\delta_{i}\log(x_{i})-\beta\sum_{i=1}^{n}x_{i}^{\alpha},$

where $r=\sum_{i=1}^{n}\delta_{i}$ is the failure number ($r<n$). Thus, making
$\dfrac{d}{d\beta}\ell(\beta,\alpha;\bm{x})=0$ and
$\dfrac{d}{d\alpha}\ell(\beta,\alpha;\bm{x})=0$, the MLE of $\beta$ is given
by:

$\hat{\beta}=\dfrac{r}{\sum_{i=1}^{n}x_{i}^{\hat{\alpha}}},$

where $\hat{\alpha}$ is the MLE of $\alpha$, which can be determine by solving
the nonlinear equation:

$\frac{r}{\alpha}+\sum_{i=1}^{n}\delta_{i}\log(x_{i})=\left(\dfrac{r}{\sum_{i=1}^{n}x_{i}^{\alpha}}\right)\sum_{i=1}^{n}x_{i}^{\alpha}\log(x_{i}).$

The MLEs for random censoring is a generalization of type I and II censoring,
therefore, they both can be achieved using the equations above. Again, the
non-linear equation has to be solved and the function uniroot is used to
obtain the estimate of $\alpha$. The codes are given by:

⬇

1likealpha<-function(alpha){

2r/alpha+sum(delta*log(t))}- (r/(sum(t^alpha)))*sum(t^alpha*log(t))

3

4###Simplifing some arguments###

5r<-sum(delta)

6

7###Computing the MLE###

8alpha<-try(uniroot(f=likealpha, c(0,100))$root)

9beta<-r/(sum(t^alpha))

The data presented in Example 4.3 is considered here, from the code above the
R output is given by

    > print(c(alpha,beta))
    [1] 1.34 2.10

The estimates for the parameters, i.e., $\hat{\alpha}=1.34$ and
$\hat{\beta}=2.10$, are not so far from the true values $\alpha=1.5$ and
$\beta=2.5$. In fact, considering $n=50,000$, the estimates for the parameter
are $\hat{\alpha}=1.50$ and $\hat{\beta}=2.51$, with $39.44\%$ of censored
data.

###  Long-term survival with random censoring

The presence of long-term survival includes an additional parameters. In this
case, a modification in the PDF and the survival function is necessary. From
the survival function (15), the long term probability density function is
given by

$f(x;\beta,\alpha,p)=\alpha\beta\left(1-p\right)x^{\alpha-1}\exp\left(-\beta
x^{\alpha}\right),$

where $\beta>0$ and $\alpha>0$ are, respectively, the scale and shape
parameters and $p\in(0,1)$ denotes the proportion of “cured” or “immune”
individuals. Let $X_{1},\ldots,X_{n}$ be a random sample of the long-term
Weibull distribution, the likelihood function considering data with random
censoring is given by:

$\displaystyle
L(\beta,\alpha,p;\bm{x,\delta})=\alpha^{r}\beta^{r}(1-p)^{r}\prod_{i=1}^{n}x_{i}^{\delta_{i}(\alpha-1)}\left[p+(1-p)\exp\left(-\beta
x_{i}^{\alpha}\right)\right]^{1-\delta_{i}}\exp\left(-\beta\sum_{i=1}^{n}\delta_{i}x_{i}^{\alpha}\right),$

where $r=\sum_{i=1}^{n}\delta_{i}$. Thus, the $\log$-likelihood function is:

$\displaystyle\ell(\beta,\alpha,p;\bm{x,\delta})=$ $\displaystyle\
r\log(\alpha)+r\log(\beta)+r\log(1-p)+(\alpha-1)\sum^{n}_{i=1}\delta_{i}\log(x_{i})-\beta\sum_{i=1}^{n}\delta_{i}x_{i}^{\alpha}$
(16)
$\displaystyle+\sum_{i=1}^{n}(1-\delta_{i})\log\left(p+(1-p)\exp\left(-\beta
x_{i}^{\alpha}\right)\right).$

The likelihood equations are given by:

$\frac{r}{\alpha}+\sum_{i=1}^{n}\delta_{i}\log(x_{i})+\eta_{1}(\beta,\alpha,p;\bm{x,\delta})=\beta\sum_{i=1}^{n}\delta_{i}x_{i}^{\alpha}\log(x_{i}),$

$\frac{r}{\beta}+\eta_{2}(\beta,\alpha,p;\bm{x,\delta})=\sum_{i=1}^{n}\delta_{i}x_{i}^{\alpha}\quad\mbox{and}\quad\eta_{3}(\beta,\alpha;\bm{x,\delta})=\frac{r}{1-p},$

where $(\theta_{1},\theta_{2},\theta_{3})=(\beta,\alpha,p)$ and

$\eta_{i}(\beta,\alpha,p;\bm{x,\delta})=\sum_{i=1}^{n}(1-\delta_{i})\dfrac{d}{d\theta_{i}}\log
S(x_{i},\beta,\alpha,p),\,\,\mbox{ for }i=1,2,3.$

The solution of the likelihood equations can be achieved following the same
approach discussed in the previous sections. On the hand, we can compute the
estimates directly from the log-likelihood equation given by (16). The codes
necessary to compute the estimates are given as follows.

⬇

1loglike<-function(theta) {

2alpha<-theta[1]

3beta<-theta[2]

4p<-theta[3]

5aux<- r*log(alpha)+r*log(1-p)+r*log(beta)+(alpha-1)*sum(delta*log(t))

6+sum((1-delta)*log(p+((1-p)*(exp(-(beta*(t^alpha)))))))

7-sum(delta*(beta*(t^alpha)))

8return(aux)

9}

10

11###Simplifing some arguments###

12r<-sum(delta)

13###Initial Values###

14alpha<-1.11

15beta<-0.88

16p<-0.6

17

18library(maxLik)

19##Calculing the MLEs

20estimate <- try(maxLik(loglike, start=c(alpha,beta,p)))

Note that we have defined a function with the log-likelihood to find the
maximum of the likelihood equation. We also have used the package maxLik
(Henningsen and Toomet, 2011). It is necessary to start the initial values to
initiate the iterative procedure. Here, only for illustrative purpose, we use
the MLEs of $\alpha$ and $\beta$ obtained from the random censoring that
returned $\tilde{\alpha}=1.11$ and $\tilde{\beta}=0.88$, while we can consider
for $p$ the proportion of censoring as initial value $\tilde{p}=0.6$. The
output of the previous code is given by

    > coef(estimate)
    [1] 1.705 4.023 0.444

From the initial values, we observe that if the presence of long-term survival
is ignored, the parameters estimates are highly biased, while the MLEs as
$\hat{\alpha}=1.71$, $\hat{\beta}=4.02$ and $\hat{p}=0.44$ with $60\%$ of
censored data. Although we have not discussed here, the maxLik function can
consider different optimization routines, constraining the parameters, and
setting the number of iterations (Henningsen and Toomet, 2011). Another
important characteristic is the easy computation of the Hessian matrix and the
standard errors used to construct confidence intervals. The matrix can be
obtained from the argument hessian, while the standard errors are obtained
directly from stdEr method.

    >print(estimate$hessian)
              [,1]       [,2]        [,3]
    [1,] -21.11378  3.2862602  -10.238921
    [2,]   3.28626 -0.8455459    3.772982
    [3,] -10.23892  3.7729819 -127.069910

    >print((-solve(estimate$hessian))^0.5)
          [,1]  [,2]  [,3]
    [1,] 0.350 0.706 0.070
    [2,] 0.706 1.841 0.246
    [3,] 0.070 0.246 0.096

    > stdEr(estimate)
    [1] 0.350 1.841 0.096

The standard errors for the estimated parameters under random censoring and
long-term survival are given above and can be used to construct the asymptotic
confidence intervals discussed in Section 2 and the coverage probabilities.

##  Simulation study

In this section, we present a Monte Carlo simulation study to illustrate
application of the pseudo-random sample schemes presented in the previous
sections. From the generate samples we estimate the parameters using the
respective MLEs. Let $\bm{\theta}$ be the parameter vector with size
$i=1,\ldots,k$ and
$\hat{\theta}_{i,1},\hat{\theta}_{i,2},\ldots,\hat{\theta}_{i,N}$ denote the
respective maximum likelihood estimates. To verify the performance of the
estimation method, the following criteria are adopted: the average bias, mean
squared error (MSE) and the coverage probabilities (CPs), which are
calculated, respectively, by

$\text{Bias}_{i}=\dfrac{1}{N}\sum_{j=1}^{N}\left(\hat{\theta}_{i,j}-\theta_{i}\right),\quad\text{MSE}_{i}=\dfrac{1}{N}\sum_{j=1}^{N}\left(\hat{\theta}_{i,j}-\theta_{i}\right)^{2}\
\ \text{and}\ \
\text{CP}(\theta_{i})=\dfrac{\sum_{j=1}^{N}I_{(a,b)}\left(\theta_{i}\right)}{N},$
(17)

where $I_{A}(x)$ means the indicator function,
$a=\hat{\theta}_{i}-z_{\frac{\xi}{2}}\sqrt{H^{-1}_{ii}(\hat{\bm{\theta}})}$
and
$b=\hat{\theta}_{i}+z_{\frac{\xi}{2}}\sqrt{H^{-1}_{ii}(\hat{\bm{\theta}})}$,
with $z_{\xi/2}$ representing the $\left(\xi/2\right)$-th quantile of a
standard normal distribution.

The Bias informs how far the average of the estimated values is from the true
parameter, whereas the MSE measures the average of the errors’ squares.
Estimators that returns both Biases and MSEs closer to zero are usually
preferred. Regarding the 95% CP, for a large number of experiments, using a
95% confidence level, the frequencies of intervals covering the true values of
$\theta_{i}$ should be closer to 95%. Although it is widely known that the
MLEs are biased for small samples, we consider this measure since the Bias
goes to zero when $n$ increases. Hence, if the algorithm to generate pseudo-
random samples with censoring does not work correctly, they should return
Biased estimates even for large samples.

###  Type-II, type-I and random censoring

We consider the MLEs discussed in Section 5 without the additional long-term
parameter. We generate N=100,000 random samples for each type-I, type-II, and
random censored scheme and its respective maximum likelihood estimates to
conduct such analysis. The steps performed by the algorithm are the following:

Algorithm 9 Weibull with type-I,type-II, and random censoring

1: Define the parameter value of $\theta=(\alpha,\beta)$;

2: Generate pseudo-random samples of size $n$ from Weibull distribution;

3: Generate the pseudo-random sample with censoring.

4: Maximize the log-likelihood function by using the MLE
$\bm{\hat{\theta}}=(\hat{\alpha},\hat{\beta})$;

5: Repeat $N$ times steps 2 till 4.

6: Compute the measures presented in (17).

The parameter values are $\alpha=1.5$ and $\beta=2.5$, while the sample sizes
are $n\in\\{10,40,60,\ldots,300\\}$ with proportion of censoring given by 0.4.
Here, we have used the same complete sample for the three different censoring
schemes, such that we can observe the effect of censoring in the Bias and MSE.
Figure 4 shows the Bias of the MSE and the coverage probability for the
parameters of the Weibull distribution.

Figure 4: MREs, MSEs and CPs related to the estimates of $\alpha=2.5$ and
$\beta=1.5$ for $N=100,000$ simulated samples, considering
$n=(10,20,\ldots,300)$.

From the graph above, we can see that as $n$ increases, both bias and MSE are
closer to zero, as expected. Moreover, the CPs are also closer to the $0.95$
nominal level. Hence, the proposed schemes to sample with censoring do not
include an additional bias than observed with the MLEs. Therefore, our
proposed approaches can be further modified for different models.

###  Long-term Weibull with random censoring

The estimation assuming the presence of a long-survival term includes an
additional parameter to be estimated. Hence, we have considered such analysis
in separate sections to account for this additional parameter. The algorithm 7
are used to obtain the pseudo-random samples with censoring, while the maximum
likelihood estimators, discussed in Section 5.4, are considered to obtain the
estimates of parameter vector $\bm{\theta}=(\alpha,\beta,p)$ from the long-
term Weibull distribution. The algorithm is described as follows:

Algorithm 10 Long-term Weibull with random censoring

1: Define the parameter value of $\theta=(\alpha,\beta,p)$;

2: Generate a random sample of size $n$ from long-term Weibull distribution;

3: Obtain the censored data using Algorithm 7;

4: Maximize the log-likelihood function by using the MLEs presented in (5.4);

5: Repeat $N$ times steps 2 till 4.

6: Compute the measures presented in (17).

Assuming the same steps as in previous subsection, we compute the average
bias, mean squared error (MSE) and coverage probability (CP) for each
$\theta_{i}$, for $i=1,2,3$. We define $n\in\\{40,50,\ldots,400\\}$ and
$N=100,000$, with $\alpha=1.5$, $\beta=2.5$ and $\pi=(0.3,0.5,0.7)$. Hence, we
expect 3 different levels of cure fraction in the population. In these cases,
the censoring proportions are, respectively $0.4$, $0.6$ and $0.8$. To
achieved such proportions we have selected $\lambda=(3.110,2.159,1.389)$,
obtained by equation (8). Figure 5 displays the obtained Bias, MSE and CP of
the estimates of the long-term survival model and random censoring data.

Figure 5: MREs, MSEs and CPs related to the estimates of $\alpha=1.5$,
$\beta=2.5$ and $\pi=c(0.3,0.5,0.7)$ for $N=100,000$ simulated samples,
considering different values of $n=(50,60,\ldots,400)$.

It can be seen from the results in Figure 5 that as $n$ increases, both bias
and MSE goes to zero, as expected. Additionally, the CPs are also converging
to the $0.95$ nominal level. The different levels of long-term survival and
censoring indicate that as the proportion of censored data increases,
uncertainty is included in the Bias, and the MSE of the estimates also
increases. On the other hand, as we have large samples, the Bias tends to the
true value regardless of the proportion of censored data.

##  Conclusions

Censored information occurs in nature, society, and engineering Klein and
Moeschberger (2006). Due to their relevance, several approaches were developed
to consider such partial information during the statistical analysis. Despite
the theoretical and applied advances in this area, the necessary mechanism to
generate pseudo-random samples has not been considered in a unified framework.
This paper has provided a unified discussion of the most common approaches to
simulate pseudo-random samples with right-censored data. The censoring schemes
considered were the most widely observed in real situations: type-I, type-II,
and random censoring. The latter censoring scheme is general in the sense that
it has the first two as particular cases.

We have presented a comprehensive discussion about the most common methods
used to obtain pseudo-random samples from probability distributions, such as
the inverse transformation method, an approach to simulate mixture models, and
the Metropolis-Hasting algorithm, which is a Markov Chain Monte Carlo method.
Further, we have provided a unified discussion on the necessary steps to
modify the pseudo-random samples to include censored information.
Additionally, we also have discussed the necessary steps to achieve a
desirable proportion of censoring. The proposed sampling schemes have been
illustrated using the Weibull distribution, one of the most common lifetime
distributions. Further, we also have considered another probability
distribution in which the sampling procedure and the estimation method are not
easily implemented, i.e., the power-law distribution with cut off. This
distribution has been considered, and the Metropolis-Hasting algorithm has
been used to obtain pseudo-random samples. The maximum likelihood estimators
have also been presented to validate the proposed approach. The codes in
language R are presented and discussed to allow the readers to reproduce the
procedures assuming other distributions.

The maximum likelihood method has been used for estimating the parameters of
the Weibull distribution by considering different types of censoring and long-
term survival. We have shown that the estimated parameters were closer to the
true values. If the samples were not generated correctly, it would be expected
that even for large samples, a systematic bias would be presented in the
estimates, which were not observed during the simulation study. The coverage
probabilities of asymptotic intervals also have presented the expected
behavior assuming that the censored information was generated correctly.

Our approach is general, and our results can be applied to most distributions
that have been introduced in the literature. These results also open new
opportunities for the analysis of censored data in models that are under
development. Although we have not discussed the presence of a regression
structure, the inverse transformed method can be used in the CDF with
regression parameters. These results can also be applied in frailty models
(Wienke, 2010) since the unobserved heterogeneity is usually included using a
probability distribution, which, combined with the baseline model, leads to a
generalized distribution. Hence, both inverse transformed methods or
Metropolis-Hasting algorithm can be used to sample complete data. Therefore,
our approach can be used to modify the data with censoring. There are a large
number of possible extensions of this current work. Other censoring
mechanisms, such as left-truncated data and the different progressive
censoring mechanisms, can be investigated using our approach; this is an
exciting topic for further research.

## Disclosure statement

No potential conflict of interest was reported by the author(s)

## Acknowledgements

Pedro L. Ramos acknowledges support from the São Paulo State Research
Foundation (FAPESP Proc. 2017/25971-0). Daniel C. F. Guzman and Alex L. Mota
acknowledge the support of the coordination for the improvement of higher-
level personnel - Brazil (CAPES) - Finance Code 001. Francisco Rodrigues
acknowledges financial support from CNPq (grant number 309266/2019-0).
Francisco Louzada is supported by the Brazilian agencies CNPq (grant number
301976/2017-1) and FAPESP (grant number 2013/07375-0).

## Apendix - MLE for the power-law with cut off

Here, we have derived the maximum likelihood estimators for the parameter of
the PLC distribution. Let $T_{1},\ldots,T_{n}$ be a random sample such that
$T$ has PDF given in (5), then the likelihood function is given by

$L(\alpha,\beta;\bm{x},x_{\min})=\frac{\beta^{n-n\alpha}}{\Gamma(1-\alpha,\beta
x_{\min})^{n}}\prod_{i=1}^{n}x_{i}^{-\alpha}\exp\left(-\beta\sum_{i=1}^{n}x_{i}\right).$
(18)

The log-likelihood function
$l(\alpha,\beta;\bm{t})=\log{L(\alpha,\beta;\bm{x},x_{\min})}$ is given by

$\displaystyle l(\alpha,\beta;\bm{x},x_{\min})=$ $\displaystyle
n(1-\alpha)\log(\beta)-\alpha\sum_{i=1}^{n}\log(x_{i})-\beta\sum_{i=1}^{n}x_{i}-n\log\Gamma\left(1-\alpha,\beta
x_{\min}\right),$ (19)

and the score functions are given by

$\displaystyle U(\alpha;\bm{x},x_{\min})=$
$\displaystyle-n\log(\beta)+\frac{G_{2,3}^{\,3,0}\\!\left(\left.{\begin{matrix}1,1\\\
0,0,1-\alpha\end{matrix}}\;\right|\,\beta
x_{\min}\right)}{\Gamma(1-\alpha,\beta x_{\min})}+\log(\beta
x_{\min})-\sum_{i=1}^{n}\log(x_{i})$ (20)
$U(\beta;\bm{x},x_{\min})=-\frac{n(1-\alpha)}{\beta}-\sum_{i=1}^{n}x_{i}+\frac{e^{-\beta
x_{\min}}}{\beta E_{\alpha}(\beta x_{\min})}.$ (21)

The maximum likelihood estimate is asymptotically normal distributed with a
bivariate normal distribution given by $(\hat{\alpha},\hat{\beta})\sim
N_{2}\left((\alpha,\beta),I^{-1}(\alpha,\beta)\right)\mbox{ for }n\to\infty$,
where $I^{-1}(\alpha,\beta)$ is the inverse of the Fisher information matrix,
where the elements are given by

$I_{\alpha,\alpha}(\alpha,\beta)=\frac{2n\Gamma(1-\alpha,\beta
x_{\min})G_{3,4}^{\,4,0}\\!\left(\left.{\begin{matrix}1,1,1\\\
0,0,0,1-\alpha\end{matrix}}\;\right|\,\beta
x_{\min}\right)}{\Gamma(1-\alpha,\beta
x_{\min})^{2}}-n\frac{G_{2,3}^{\,3,0}\\!\left(\left.{\begin{matrix}1,1\\\
0,0,1-\alpha\end{matrix}}\;\right|\,\beta
x_{\min}\right)^{2}}{\Gamma(1-\alpha,\beta x_{\min})^{2}},$
$I_{\alpha,\beta}(\beta,\alpha)=I_{\beta,\alpha,}(\alpha,\beta)=\frac{n}{\beta}-nx_{\min}e^{-\beta
x_{\min}}\left[\frac{G_{2,3}^{\,3,0}\\!\left(\left.{\begin{matrix}1-\alpha,1-\alpha\\\
1-2\alpha,-\alpha,-\alpha\end{matrix}}\;\right|\,\beta
x_{\min}\right)}{\Gamma(1-\alpha,\beta x_{\min})^{2}}\right],$ $\displaystyle
I_{\beta,\beta}(\alpha,\beta)=-\frac{ne^{-2\beta x_{\min}}\left(1-e^{\beta
x_{\min}}\left(\alpha+\beta x_{\min}\right)E_{\alpha}\left(\beta
x_{\min}\right)\right)}{\beta^{2}E_{\alpha}\left(\beta
x_{\min}\right)^{2}}-\frac{n\left(\alpha-1\right)}{\beta^{2}}.$

The score function and the elements of the Fisher information matrix depend on
Meijer-G functions which are not easy to be computed. On the other hand, we
can obtain the estimates directly from the maximization of the loglikelihood
function. Additionally, we observe that observed Fisher information is the
same as the expected Fisher information matrix since they do not depend on
$x$, hence we can construct the confidence intervals directly from the Hessian
matrix that is computed from the maxLik function. The necessary codes to
compute the results above are displayed in:

⬇

1#The incomplete gamma function

2ig <- function(x,a) {x^(a-1)*exp(-x)}

3#The loglikelihood function

4loglike <- function(theta){

5 alpha <- theta[1]

6 beta <- theta[2]

7 l<- n*(1-alpha)*log(beta)-alpha*sum(log(x))-beta*sum(x)

8 -n*log(integrate(ig, lower = (beta*xm), upper = Inf,a=(1-alpha))$value)

9 return(l) }

10

11#Computing the MLEs and the Hessian matrix

12res <- maxLik(loglike, start=c(0.5, 0.5))

As initial values, we have used $\tilde{\alpha}=0.5$ and $\tilde{\beta}=0.5$.
The output containing the parameter estimates as well as its related variances
are given by:

    > coef(res)
    [1] 1.655 0.281
    > stdEr(res)
    [1] 0.428 0.140

From the results above, we have that $\hat{\alpha}=1.655$ and
$\hat{\beta}=0.281$ with the respective standard errors $0.428$ and $0.140$.
The confidence intervals is construct from the asymptotic theory and are
$(0.817;2.492)$ for $\alpha$ and $(0.007;0.555)$ for $\beta$.

## References

  * Algarni et al. (2020) Algarni, A., A. M. Almarashi, H. Okasha, and H. K. T. Ng (2020).  E-bayesian estimation of chen distribution based on type-i censoring scheme.  Entropy 22(6), 636.
  * Balakrishnan et al. (2007) Balakrishnan, N., D. Kundu, K. T. Ng, and N. Kannan (2007).  Point and interval estimation for a simple step-stress model with type-ii censoring.  Journal of Quality Technology 39(1), 35–47.
  * Berkson and Gage (1952) Berkson, J. and R. P. Gage (1952).  Survival curve for cancer patients following treatment.  Journal of the American Statistical Association 47(259), 501–515.
  * Bogaerts et al. (2017) Bogaerts, K., A. Komarek, and E. Lesaffre (2017).  Survival analysis with interval-censored data: A practical approach with examples in R, SAS, and BUGS.  CRC Press.
  * Chen et al. (1999) Chen, M.-H., J. G. Ibrahim, and D. Sinha (1999).  A new bayesian model for survival data with a surviving fraction.  Journal of the American Statistical Association 94(447), 909–919.
  * Chen et al. (2002) Chen, M.-H., J. G. Ibrahim, and D. Sinha (2002).  Bayesian inference for multivariate survival data with a cure fraction.  Journal of Multivariate Analysis 80(1), 101–126.
  * Chib and Greenberg (1995) Chib, S. and E. Greenberg (1995).  Understanding the metropolis-hastings algorithm.  The american statistician 49(4), 327–335.
  * Clauset et al. (2009) Clauset, A., C. R. Shalizi, and M. E. Newman (2009).  Power-law distributions in empirical data.  SIAM review 51(4), 661–703.
  * Dellaportas et al. (1997) Dellaportas, P., D. Karlis, and E. Xekalaki (1997).  Bayesian analysis of finite poisson mixtures.  Athens University of Economics and Business, Greece.
  * Feller (2008) Feller, W. (2008).  An introduction to probability theory and its applications, Volume 2.  John Wiley & Sons.
  * Field et al. (2012) Field, A. P., J. Miles, and Z. Field (2012).  Discovering statistics using r.
  * Gamerman and Lopes (2006) Gamerman, D. and H. F. Lopes (2006).  Markov chain Monte Carlo: stochastic simulation for Bayesian inference.  Chapman and Hall/CRC.
  * Gelman et al. (1992) Gelman, A., D. B. Rubin, et al. (1992).  Inference from iterative simulation using multiple sequences.  Statistical science 7(4), 457–472.
  * Ghosh and Basu (2019) Ghosh, A. and A. Basu (2019).  Robust and efficient estimation in the parametric proportional hazards model under random censoring.  Statistics in Medicine 38(27), 5283–5299.
  * Guzmán et al. (2020) Guzmán, D. C., C. S. Ferreira, and C. B. Zeller (2020).  Linear censored regression models with skew scale mixtures of normal distributions.  Journal of Applied Statistics, 1–26.
  * Hastings (1970) Hastings, W. K. (1970).  Monte carlo sampling methods using markov chains and their applications.
  * Henningsen and Toomet (2011) Henningsen, A. and O. Toomet (2011).  maxlik: A package for maximum likelihood estimation in r.  Computational Statistics 26(3), 443–458.
  * Hoel et al. (1986) Hoel, P. G., S. C. Port, and C. J. Stone (1986).  Introduction to stochastic processes.  Waveland Press.
  * James et al. (2013) James, G., D. Witten, T. Hastie, and R. Tibshirani (2013).  An introduction to statistical learning, Volume 112.  Springer.
  * Kim et al. (2011) Kim, C., J. Jung, and Y. Chung (2011).  Bayesian estimation for the exponentiated weibull model under type ii progressive censoring.  Statistical Papers 52(1), 53–70.
  * Kim and Jhun (2008) Kim, Y.-J. and M. Jhun (2008).  Cure rate model with interval censored data.  Statistics in medicine 27(1), 3–14.
  * Klein and Moeschberger (2006) Klein, J. P. and M. L. Moeschberger (2006).  Survival analysis: techniques for censored and truncated data.  Springer Science & Business Media.
  * Li (2003) Li, L. (2003).  Non-linear wavelet-based density estimators under random censorship.  Journal of statistical planning and inference 117(1), 35–58.
  * Maller and Zhou (1996) Maller, R. A. and X. Zhou (1996).  Survival analysis with long-term survivors.  John Wiley & Sons.
  * Martinez et al. (2016) Martinez, E. Z., J. A. Achcar, M. V. de Oliveira Peres, and J. A. M. de Queiroz (2016).  A brief note on the simulation of survival data with a desired percentage of right-censored datas.  Journal of Data Science 14(4), 701–712.
  * Mehta et al. (2019) Mehta, P., M. Bukov, C.-H. Wang, A. G. Day, C. Richardson, C. K. Fisher, and D. J. Schwab (2019).  A high-bias, low-variance introduction to machine learning for physicists.  Physics Reports 810, 1–124.
  * Metropolis et al. (1953) Metropolis, N., A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller (1953).  Equation of state calculations by fast computing machines.  The journal of chemical physics 21(6), 1087–1092.
  * Migon et al. (2014) Migon, H. S., D. Gamerman, and F. Louzada (2014).  Statistical inference: an integrated approach.  CRC press.
  * Niederreiter (1992) Niederreiter, H. (1992).  Random number generation and quasi-Monte Carlo methods.  SIAM.
  * Peng and Dear (2000) Peng, Y. and K. B. Dear (2000).  A nonparametric mixture model for cure rate estimation.  Biometrics 56(1), 237–243.
  * Piegorsch (1990) Piegorsch, W. W. (1990).  Maximum likelihood estimation for the negative binomial dispersion parameter.  Biometrics, 863–867.
  * Qin and Jing (2001) Qin, G. and B.-Y. Jing (2001).  Empirical likelihood for cox regression model under random censorship.  Communications in Statistics-Simulation and Computation 30(1), 79–90.
  * R Core Team (2019) R Core Team (2019).  R: A Language and Environment for Statistical Computing.  Vienna, Austria: R Foundation for Statistical Computing.
  * Robert et al. (2010) Robert, C. P., G. Casella, and G. Casella (2010).  Introducing monte carlo methods with r, Volume 18.  Springer.
  * Rodrigues et al. (2009) Rodrigues, J., V. G. Cancho, M. de Castro, and F. Louzada-Neto (2009).  On the unification of long-term survival models.  Statistics & Probability Letters 79(6), 753–759.
  * Ross et al. (2006) Ross, S. M. et al. (2006).  A first course in probability, Volume 7.  Pearson Prentice Hall Upper Saddle River, NJ.
  * Schumacker and Tomek (2013) Schumacker, R. and S. Tomek (2013).  Understanding statistics using R.  Springer Science & Business Media.
  * Shilo et al. (2020) Shilo, S., H. Rossman, and E. Segal (2020).  Axes of a revolution: challenges and promises of big data in healthcare.  Nature Medicine 26(1), 29–38.
  * Singh et al. (2005) Singh, U., P. K. Gupta, and S. K. Upadhyay (2005).  Estimation of parameters for exponentiated-weibull family under type-ii censoring scheme.  Computational statistics & data analysis 48(3), 509–523.
  * Suess and Trumbo (2010) Suess, E. A. and B. E. Trumbo (2010).  Introduction to probability simulation and Gibbs sampling with R.  Springer Science & Business Media.
  * Tan et al. (2016) Tan, P.-N., M. Steinbach, and V. Kumar (2016).  Introduction to data mining.  Pearson Education India.
  * Tse et al. (2000) Tse, S. K., C. Yang, and H.-K. Yuen (2000).  Statistical analysis of weibull distributed lifetime data under type ii progressive censoring with binomial removals.  Journal of Applied Statistics 27(8), 1033–1043.
  * Tsodikov et al. (1996) Tsodikov, A. D., A. Y. Yakovlev, and B. Asselain (1996).  Stochastic models of tumor latency and their biostatistical applications, Volume 1.  World Scientific.
  * Wang and Wang (2001) Wang, Q. and J.-L. Wang (2001).  Inference for the mean difference in the two-sample random censorship model.  Journal of Multivariate Analysis 79(2), 295–315.
  * Wang and Li (2002) Wang, Q.-H. and G. Li (2002).  Empirical likelihood semiparametric regression analysis under random censorship.  Journal of Multivariate Analysis 83(2), 469–486.
  * Weibull (1939) Weibull, W. (1939).  A statistical theory of strength of materials.  IVB-Handl..
  * Wienke (2010) Wienke, A. (2010).  Frailty models in survival analysis.  CRC press.
  * Yin and Ibrahim (2005) Yin, G. and J. G. Ibrahim (2005).  Cure rate models: a unified approach.  Canadian Journal of Statistics 33(4), 559–570.