% !TEX TS-program = LuaLaTeX
% !TEX encoding = UTF-8 Unicode 
% !TEX root = ../../mal-core.tex

\def\d{\mbox{\rm d}}
\def\KM{\textrm{KM}}

\renewcommand{\arraystretch}{1.3}
\renewcommand\citeform[1]{}
\def\citeleft{}
\def\citeright{}


\gdef\mujnazevCS{Přehled vybraných metod řešení úlohy\odskoc lineární regrese s~cenzorovanými daty}
%P\v{R}EHLED VYBRAN\'YCH METOD \v{R}E\v{S}EN\'I \'ULOHY\\
%LINE\'ARN\'I REGRESE S CENZOROVAN\'YMI DATY}
\gdef\mujnazevEN{An overview of methods of linear regression analysis with censored data}
\gdef\mujnazevPR{\mujnazevCS}
\gdef\mujnazevDR{\mujnazevEN}
\gdef\mujauthor{Petr Volf}

\bonushyper

\nazev{Přehled vybraných metod řešení úlohy\\[2pt] lineární regrese s~cenzorovanými daty}

\nazev{An overview of methods of linear regression\\[2pt] analysis with censored data}

\author{\mujauthor}

\Adresa{ Faculty of Sciences, Humanities and Education,
   Technical University of Liberec, Czech Republic}
   
\Email{petr.volf@tul.cz}


\Abstrakt{Po úspěšném rozšíření a~používání Coxova regresního modelu
v~70. letech (a~někdy až jeho nadužívání bez ohledu na splnění předpokladu proporcionality rizik) se pochopitelně zvýšila snaha najít vhodné metody také pro analýzu lineárního regresního modelu s~cenzorovanými daty. Tento příspěvek přibližuje tři základní ``historické'' metody (autorů Miller, 1976; Buckley and James, 1979; Koul et al., 1981) a~vzájemně je porovnává (na jednoduchém příkladu s~náhodně generovanými daty).
Výsledky jsou porovnány i~s~řešením pomocí maxima věrohodnosti 
při předpokladu normálně rozdělených odchylek, i~s~obecnějším 
semi-parametrickým modelem se zrychleným časem (AFT model).}
\smallskip
\KlicovaSlova{Lineární regresní model, metoda nejmenších čtverců, cenzorovaná data, AFT model, Coxův model}

\Abstract{As a~reaction to results of D.\,R. Cox and others in '70s
proposing regression models able to work with censored data,
a~number of statisticians (e.g. Miller, 1976; %\cite{Mill}%
Buckley and James, 1979; %\cite{BuJa}%
Koul et al., 1981) %\cite{KSVR}%
have tried to extend the method of the
least-squares in the same sense, i.e. to adapt it to censoring.
The present contribution contains an overview of these approaches,
compares them with direct numerical solution of corresponding
maximal likelihood task and recalls also the AFT (Accelerated
Failure Time) model as a~semi-parametric generalization of the
linear regression model. Performance of methods will be shown on
simple numerical examples.}
\smallskip
\KeyWords{Linear regression model, method of least squares, censoring,
AFT model, Cox model}



\section{Introduction}

The standard linear regression model for $N$ pairs of variables
$(X_i, Y_i), \,i = 1, ...,N$ has the form
 $$
Y_i = \alpha +\beta \cdot X_i + \varepsilon_i,\,\quad i=1,...,N,
\eqno{(1)}
 $$
with $\varepsilon_i$ assumed to be the i.i.d. variables with zero
mean. The primary type of censorship considered here will be the
random censoring from the right side. It means that there exist
random variables $C_i$, say, independent mutually and of $Y_i$,
and instead $Y_i$ the values of $Y^*_i=\min (Y_i,C_i)$ are
observed. As a~rule,  it is known whether an observation is censored
or not, which is denoted with an indicator
$\delta_i=1$ if $Y_i\le C_i,\,\, \delta_i=0$ if $Y_i>C_i$.

After the introduction of the Cox model (Cox, 1972), 
\cite{Cox}
the effort of
statisticians has concentrated also to the solution of linear
regression problem with censoring. In the present overview we
shall recall several most important and more or less successful
approaches provided in Miller (1976), \cite{Mill}%
Buckley and James (1979), \cite{BuJa}%
and Koul et al. (1981). \cite{KSVR}%
It is also interesting to follow these
methods development from the point of view of their assumptions
and of their real performance. Let us summarize first certain
common assumptions:
\begin{itemize}
\itemsep=0pt
\item[1.] Residuals $\varepsilon_i$ are assumed to be i.i.d. with
zero mean and finite variance.

\item[2.] Censoring variables $C_i$ are independent
mutually and of $\varepsilon_i$.

\item[3.] In general, $C_i$ can depend on covariates $X_i$, however
some approaches do not allow it, namely that of Koul et al. 
\cite{KSVR}%
There
the i.i.d. $C_i$-s are required.

\item[4.] When the estimation is based on the maximum likelihood
principle, it has to be assumed that variables $C_i$ do not
contain any information on model parameters (this is called
``non-informative'' censoring, see also the assumptions for
estimation in the Cox regression model); only then we can omit the
likelihood part containing distribution of $C_i$.
\end{itemize}

This contribution recapitulates the above mentioned
approaches first. Then the approach based on the MLE is presented
(based on assumed type of the distribution of residuals).
Further, the relationship of the linear model with
the Accelerated Failure Time (AFT) regression model is shown.
The performance of all methods is
illustrated and compared with the aid of an example. Finally,
also the Cox regression model is recalled, in order to show
a~connection of its parametrized version (with Weibull baseline
distribution) with the linear regression problem

\section{Overview of methods}

All three traditional methods described bellow have one common
feature, they all use the Kaplan-Meier ``Product Limit Estimate''
(PLE), though in slightly different ways. As the PLE is
a~nonparametric estimator of distribution function (from censored
data), its use may lead to a~lower precision of estimation (in
comparison to the case without censoring). Simultaneously, the
analysis of asymptotic properties of estimates is then more
difficult, confidence intervals for parameters can be obtained, as
a~rule, just with the aid of bootstrap.

Let us recall here briefly several basic notions from the area of
statistical survival analysis: Consider a~positive continuous-type
random variable~$Z$, with density $f$ and distribution function
$F$. $\overline{F}=1-F$ is then called the survival function, the
hazard rate is defined as $h(z)=-\d \ln (\overline{F}(z))/\d z =
f(z)/\overline{F}(z)$, the cumulative hazard rate (CHR) then equals
$H(z)=\int_0^z h(t)\d t = -\ln (\overline{F}(z))$.

The form of the Kaplan-Meier PLE of distribution function $F$ can
be found elsewhere (for instance in the most of references cited
here): Let $Z^*_i$ be ordered sample of randomly right-censored
realizations of $Z$, i.e. $Z^*_1\le Z^*_2 \le ... \le Z^*_N,$
$\delta_i$ corresponding censoring indicators, then the PLE equals
 $$
 F^{\KM}(z)=1 - \prod_{Z^*_i\le z} {\Big (}{N-i\over
 N-i+1}{\Big)}^{\delta_i}.
 $$
 Let us mention here also one practical detail
of the PLE construction: the maximal value of data (here $Z^*_N$)
is always taken as non-censored.

\subsection{Miller's estimator}

The way of construction of the Miller's estimator (Miller, 1976)
\cite{Mill}%
seems to be rather reasonable, though there also is a~common
problem with proving estimator properties (its consistency and
asymptotic normality). In particular, Miller has just shown his
estimator being consistent when the censoring variable fulfils the
regression with the same slope, i.e. $C\sim \alpha^c+\beta\cdot
X$.

The usual least squares estimates $a,b$ of $\alpha, \beta$ of
model (1) are those values which minimize the sum of squares
$\sum_{i=1}^N \varepsilon_i(a,b)^2=\sum_{i=1}^N (Y_i - a - b
X_i)^2.$ This is equivalent to minimizing $\sum
\varepsilon_i(a,b)^2\cdot\Delta \widehat{F}(\varepsilon_i(a,b))$,
where $\widehat{F}$ is the empirical distribution function (EDF) of
residuals $\varepsilon_i(a,b)=Y_i-a- bX_i$ for fixed $a,b$,
$\Delta \widehat{F}$ are the increments of $\widehat{F}$. Residuals
$\varepsilon_i(a,b)$ are available just partially, with the same
censoring pattern as $Y^*_i$, i.e. for given $a,b$ we can compute
$\varepsilon^*_i(a,b)=Y^*_i-a- bX_i$ and the censoring indicators
are the same $\delta_i$. Therefore, a~natural step is to use the
PLE of distribution function of residuals generalizing its
standard EDF to the censoring case. Namely, the task is to
minimize, over $a,b$,
 $$
\sum \varepsilon^*_i(a,b)^2\cdot\Delta
F^{\KM}(\varepsilon_i^*(a,b)), \eqno{(2)}
 $$
where $F^{\KM}$ is now the PLE of the distribution of residuals
$\varepsilon(a,b)$ for given $a,b$ (the increments $\Delta F^{\KM}$
are 0 at censored values).

There is no problem to find optimal solution numerically, with the
aid of a~convenient search method. Further, notice that for fixed
$b$ the minimum with respect to $a$ can be computed directly
without any iteration, $a$ just shifts the PLE. Similarly like in
the standard case, the minimizing $a$ equals the average of
$Y^*_i-bX_i$ w.r. to the PLE, namely
 $$
\widehat{a_b}=\sum \varepsilon^*_i(0,b)^2\cdot\Delta
F^{\KM}(\varepsilon_i^*(0,b)), \eqno{(3)}
 $$
with $\varepsilon_i^*(0,b)=Y^*_i-bX_i$. In such a~way the search
can be divided into two iterated steps, each solving just an
one-dimensional problem.

\smallskip
Unlike Miller's, the following two methods are based on attempts
to reconstruct censored data first, then these reconstructed
values are used in the ordinary least squares.

\subsection{Buckley and James}

This method of analysis in the linear regression model with
right-censoring, in fact the most popular one, has been proposed
by Buckley and James (1979). 
\cite{BuJa}%
It consists in a~version of the EM
(Expectation-Maximization) algorithm. Namely, the E step
reconstructs censored values via conditional expectation, the M
step then recomputes the model parameters from these values, using
the least squares method.

Let the model have again the form (1). Again, no specific
parametric form is assumed for the distribution of the error terms
$\varepsilon_i$. They are just assumed to be i.i.d. Buckley and
James have proposed to use weighted observations
 $$
Y^w_i = Y_i\cdot \delta_i + E(Y_i| Y_i >C_i)\cdot (1-\delta_i).
\eqno{(4)}
 $$
It is seen that uncensored values are preserved, while censored
are ``corrected'', increased. It is easy to show that for each $x$
$E(Y^w|X=x) = E(Y|X=x)$:

 \smallskip
 For each  fixed $x$ (let us omit it for the moment) and for
 each fixed $C$ we have that
 $$
 E(Y^w|C)=E(Y\cdot\delta|C)+E[E(Y|Y>C)\cdot(1-\delta)|C] =
 $$
 $$
  =\int_{-\infty}^C {y\cdot \d F_y(y)\over F_y(C)}\cdot P(\delta=1|C)
  +\int_C^{\infty} {y\cdot \d F_y(y)\over 1-F_y(C)}\cdot P(\delta=0|C)  =
 $$
 $$
 =\int_{-\infty}^C y\cdot \d F_y(y)+\int_C^{\infty} y\cdot \d F_y(y) = EY,
 $$
independently of $C$, hence $E(Y^w)=E(Y)$ as well. Here $F_y$
denotes the distribution function of $Y$ for given fixed $X=x$.

The procedure of estimation starts from an initial estimate of
$a^{(1)},\,b^{(1)}$, obtained for instance by the least squares
method from values $(Y^*_i, X_i)$, and computes corresponding
residual values $\varepsilon^*_i$, which again represent censored
values of $\varepsilon_i(a^{(1)},b^{(1)})$,  with preserved
indicators $\delta_i$. The nonparametric estimator of their
distribution, $F^{\KM}$, is again obtained via the Kaplan-Meier
PLE. From it we can ``reconstruct'' censored residuals,
 $$
\varepsilon^{(1)}_i={E}(\varepsilon_i
|\varepsilon>\varepsilon^*_i)=
{\sum_{\varepsilon^*_j>\varepsilon^*_i} \Delta
F^{\KM}(\varepsilon^*_j) \cdot \varepsilon^*_j \over
1-F^{\KM}(\varepsilon^*_i)}, \eqno{(5)}
 $$
Hence, censored outputs are then reconstructed as well, namely
 $$
Y^{(1)}_i=Y_i\cdot \delta_i + (a^{(1)}+b^{(1)} \cdot X_i +
\varepsilon^{(1)}_i)(1-\delta_i). \eqno{(6)}
 $$
Then new estimates, $a^{(2)}, b^{(2)}$ are computed via the least
squares from innovated outputs $Y^{(1)}_i$. Quite similarly,
 $$
Y^2_i=Y_i\cdot \delta_i + (a^{(2)}+b^{(2)} \cdot X_i +
\varepsilon^{(2)}_i)(1-\delta_i),
 $$
where $\varepsilon^{(2)}_i = Y_i -a^{(2)}-b^{(2)} \cdot X_i$.
 This procedure could be repeated several times till convergence.

Later on, James and Smith (1984) 
\cite{JaSm}
have established the conditions
and arguments for the consistency of the BJ estimator. However,
the problem lies in the practical computation, in its iterative
scheme. The convergence is not guaranteed, it is known that the
iterations may end in oscillation between two or more set of
values (see for instance Jin et al., 2006).
\cite{Jin}%

\subsection{Koul, Susarla, and Van Ryzin}

While the Buckley--James estimator uses (repeatedly) reconstructed
censored values, the estimator of Koul et al. (1981) 
\cite{KSVR}
reconstructs
all values, creates ``synthetic data''. This does not seem to be
an ideal concept, as it changes the basic features of the data.
Namely, they propose, instead of original censored data, to use
 $$
Y^K={\delta\cdot Y^*\over \overline{G}(Y^*)}, \eqno{(7)}
 $$
where $\overline{G}(.)$ denotes the survival function of the censoring
variable $C$, which, if unknown, is substituted by its
Kaplan-Meier PL estimator. A~consequence of such a~transformation
is that the regression of $Y^K$ on $X$ is strongly
heteroscedastic; notice that all censored values are now set to
zero, while uncensored values are increased, sometimes rather
considerably. Nevertheless, it still holds that
$E(Y^K|X=x)=E(Y|X=x)$. It could again be proved easily:

 \smallskip
In the model (1) denote $F()$ the
distribution function of errors $\varepsilon$. Then
 $$
 E(\delta\cdot Y^*)=\int_{-\infty}^{\infty} t \overline{G}(t)
 \, \d F(t-\alpha-\beta x), \mbox{    hence}
 $$
 $$
E\left({\delta\cdot Y^*\over \overline{G}(Y^*)}\right)=\int_{-\infty}^{\infty} t
 \,\d F(t-\alpha-\beta x)=\alpha+\beta x.
$$
However, as the estimation method is based on the least squares
(using reconstructed data), in fact its weighted variant should be
used. Even the authors have pointed out that the variability of
synthetic data around the response curve is neither constant nor
symmetric and may expand, therefore the estimation of response
curve may not be reliable. Figure 3 of Examples shows a~result of
such a~``data reconstruction''.

\section{A~general formulation of likelihood with censored data}

All approaches described above utilize the method of least
squares. This method corresponds in fact to the maximum likelihood
estimation under the assumption of normally distributed errors.
Therefore, the solution based on normal distribution of errors
should be equivalent, moreover it avoids the use of the
non-parametric PLE.

On the other hand, in general, we need not to restrict ourselves
just to normal model, there are, especially in the context of
reliability data studies, other natural possibilities. The linear
model can arise from the logarithmization of typical
``reliability'' models for a~variable $T$ (characterizing for
instance the time to failure) described for instance by the
Weibull, log-normal, log-logistic, gamma, or other similar
distributions. Namely, let again a~variable $Y$ fulfil the linear
regression model, $Y=\beta\cdot X + \varepsilon$. Denote
$T=\exp(Y)$ and $T_0=\exp(\varepsilon)$, their distribution
functions by $F_T,\, F_0$, resp. Then
 $$
 F_T(t)=P(T\le t)= P(\exp(\beta\cdot X)\cdot T_0 \le t)
 =F_0(t\cdot \exp(-\beta\cdot X)),
  \eqno{(8)}
 $$
which corresponds to the parametric form of the AFT model. And
vice versa, natural choices for distribution of $\varepsilon$ in
the linear regression model are then the Gumbel
(doubly-exponential, obtained by the logarithmization of the
Weibull distribution), normal, logistic, exp-gamma, and other
similar distributions.

Thus, depending on assumed distribution of errors $\varepsilon_i$,
a~direct maximization of the log-likelihood function corresponding
to this distribution, with the aid of a~convenient optimization
procedure, could be a~good choice. Let us recall that the
``relevant'' log-likelihood in the right-censored data case has the
form (provided the censoring is ``noninformative'', only then the
part concerning the distribution of variables $C$ may be omitted):
 $$
 L=\sum_{i=1}^N [\delta_i\cdot \ln(f_i(Y^*_i,\theta))+(1-\delta_i)\cdot
  \ln(1-F_i(Y^*_i,\theta))], \eqno{(9)}
  $$
where $f_i\,F_i$ now denote the density and distribution function
of variables $Y_i$, hence they depend on covariates, $\theta$ here
denotes the model parameters. The MLE, via solving equation
$\d L/\d\theta =0$, can sometimes be found quite easily (in several
steps of the Newton-Raphson algorithm); this is, among others, the
instance of the Gumbel (log-Weibull) distribution, while the case
with normal distribution has to be solved with the aid of a~more
clever iterative search procedure including repeated evaluation of
normal distribution function.

\subsection{Semi-parametric AFT regression model}

 Let us recall here also a~general version of this model.
For more information see e.g. Bagdonavi\v{c}ius and Nikulin (2002).
\cite{BaNi,KaPr}%
The accelerated failure time model is often considered as an
alternative to the popular proportional hazard model of Cox, when
the proportionality of hazards does not hold. The model assumes
that the individual speed of time is changed by a~factor depending
on covariates. Standardly this factor has the form $\exp(\beta
x)$, where $x$ is a~covariate, mostly constant in time. It follows
that the distribution of time to failure $T_i$ of an object with
covariate value $x$ has the distribution function $F(t)=F_0(t\cdot
\exp(\beta x))$, where $F_0$ characterizes a~baseline distribution
(of a~random variable $T_0$ with covariate $x=0$). It also means
that $T_{0i}=T_i\cdot \exp(\beta x_i)$ is an i.i.d. representation
of $T_0$. The logarithmic transformation of the model then leads
to the linear regression model as shown in (8). Notice also
changed sign of parameter $\beta$ in both models.

Hence, the statistical inference in the AFT model setting has to
deal with unknown baseline distribution of $T_0$, the analysis is
often complicated by the presence of censored data. As we have
seen, the estimation in a~linear model with censoring is not quite
easy task (especially when the distribution of errors is not
known). Therefore, similarly like in the case of Cox model, the
approach based on hazard rates could be preferred. Namely, let
$\{T^*_i, x_i, \delta_i,$ $i=1,..,N\}$ be observed times of
failures or censoring of $N$ objects, their covariates, indicators
of censoring, respectively. Then the likelihood equals (compare it
with (9)):
 $$
 L = \prod_{i=1}^N  h_i(T^*_i)^{\delta_i} \cdot \exp
 \left(-\int_0^{T^*_i} h_i(t)
 \d t\right) , \eqno{(10)}
 $$
where $h_i(t)=h_0(t\cdot \exp(\beta x_i))\cdot \exp(\beta x_i)$ is
the hazard rate of $i$-th object at time $t$, $h_0$ is the
baseline hazard rate of $T_0$. Theory of estimation and asymptotic
properties are derived in Lin et al. (1998) 
\cite{LiWe}
and further developed
also in Bagdonavi\v{c}ius and Nikulin (2002, Ch. 6). 
\cite{BaNi}%
A~Bayes approach
offers a~reasonable alternative, too (see Volf and Timkov\'a,
2014).
\cite{VoTi}%

Lin et al. (1998) 
\cite{LiWe}%
have shown that instead of the exact (and rather
complicated) score function for $\beta$ (i.e. obtained by the
derivation of the log-likelihood w.r. to $\beta$), it is possible
to use approximate score functions. Namely, the score function has
the form
 $$
 U(\beta)= \sum_{i=1}^N  \left\langle x_i W_i(s) -
{ \sum_j x_j W_j(s) I_j^*(T_{0i})\over \sum_k I_k^*(T_{0i})}
\right\rangle \delta_i,
 \eqno{(11)}
  $$
where $I_k^*(t)$ are indicators of risk in the scale of
$T_{0i}=T_i\cdot \exp(\beta x_i)$, hence the unknown parameter
$\beta$ is hidden in them. While exact weights $W_i(s)$ depend on
$h_0(s)$ and also on its first derivative, they may be substituted
by a~set of simpler functions, among them also by $W_i(t)=1$ for
all $i$ and~$t$. Lin et al. (1998) 
\cite{LiWe}%
have proved that corresponding
estimator of $\beta$ retains good asymptotic properties. Hence,
from such a~score function it is possible to estimate $\beta$
without knowledge of $h_0(t)$, similarly like with the aid of the
partial likelihood in the Cox model case. More details can be
found also in Nov\'ak (2013), 
\cite{Nov}%
who has proposed a~method of
goodness-of-fit test based on the random generation of residual
processes, and has studied the test behavior in various
situations.


\begin{figure}[!b]
\centering
\includegraphics[width=\textwidth]{Volf_f1.eps}
\caption{Values of $Y$ before censoring (left) and after, censored
values of $Y^*$ are denoted by `o' (right).}%
\label{volf1}%
\end{figure}

\begin{figure}[!ht]
\centering
\includegraphics[width=\textwidth]{Volf_f2.eps}
\caption{Results of the BJ procedure after 2nd iteration.
Reconstructed censored values are denoted by `o'.}%
\label{volf2}%
\end{figure}




\section{Examples}

We shall use here the following configuration of input variables
and parameters: In the LRM (1), let us set the number of data
$N=100$, the covariates $X_i$ distributed randomly uniformly in
(0,2), regression parameter $\beta=-1$, intercept $\alpha=2$,
residuals $\varepsilon_i$ following the normal distribution with
$\mu=0,\,\sigma=0.1$. Further, let the censoring values be selected
uniformly in the interval $(0,3)$, the rate of censoring is then about
30\%.

\subsection{Example 1}

First, just one set of data was generated, model parameters were estimated
from them with the aid of all methods described above.
Figure~\ref{volf1} shows the data before and after censoring.

From the non-censored data, estimates were
$a= 2.0282\,(1.9911,2.0653)$,
$b=-1.0269\, (-1.0596,-0.9942)$, with 95\% confidence intervals in
parentheses, standard deviation of residuals was estimated as
$s_r=0.0978$. Optimal line is displayed in Figure 1, left side.
For comparison, parameters were estimated also from censored data,
not taking censoring into account, in order to show the bias of
such a~`naive' approach. The following estimates were obtained:
$a=1.6172,\,b=-0.7868,\,s_r=0.3568$, corresponding line is shown
in the right part of Figure~\ref{volf1}. These values were later used as
initial estimates in the BJ procedure.

\begin{figure}[p]
\centering
\includegraphics[width=.7\textwidth]{Volf_f3.eps}
\caption{Values $Y^K(X)$ reconstructed via the Koul et al. method
(covariate $x$ is on horizontal axis). All censored values,
denoted by `o', were set to zero.}%
\label{volf3}%
\end{figure}

\begin{figure}[p]
\centering
\includegraphics[width=.7\textwidth]{Volf_f4.eps}
\caption{AFT model: Nonparametric estimate of the baseline CHR
(full curve) and its comparison with `true' log-normal CHR (dashed
curve).}%
\label{volf4}%
\end{figure}


Then, the methods of estimation described above were tested on the
same data, with the following results:

\begin{itemize}
\itemsep=-1pt
\item[1.] {\it Miller's estimator} yielded
$a=1.904,\,b=-0.936,\,s_r=0.165$.

\item[2.] {\it Buckley and James method:}
Figure~\ref{volf2} presents the result of 2nd iteration of their procedure.
It was revealed that there the smallest standard deviation of
residuals was achieved, $s_r=0.0912$, corresponding estimates of
parameters were $a=2.0645,\,b= -1,0415.$ After that, the results
of further steps started to oscillate.

\item[3.] {\it Koul et al. estimator:} Figure~\ref{volf3} shows the data after
proposed way of ``reconstruction". It illustrates rather clearly
the drawback of this method. Nevertheless, in this instance, the
estimates were quite reasonable, $a=1.9228,\,b=-0.9696$. However,
as expected, estimated $s_r=0.3731$ was much higher.

\item[4.] {\it Direct MLE} under the assumption of normally distributed
residuals yielded $a=2.0208,\,b=-1.0103,\,s_r=0.9853$.

\item[5.] {\it Semi-parametric estimator} in the framework of the
AFT model: We are able to estimate just the slope parameter
$\beta$, while the rest of model is hidden in the non-parametric
estimate of the cumulative baseline hazard rate. In our case, this
estimate should be close to the CHR of the log-normal distribution
with parameters $\mu=\alpha=2,\,\sigma=0.1$. The proximity of both curves
is shown in Figure~\ref{volf4}. Estimate of $\beta$ was $b=-1.0369$.
\end{itemize}




\subsection{Example 2}

The purpose was to explore the reliability of all methods in more
detail. In particular, we were interested in the precision of
estimates of the regression parameter $\beta$. Therefore, each
simulation (in the framework of the same model configuration as in
Example 1) was repeated 200 times, for each method. In such a~way,
from each method, a~set of 200 estimates of $\beta$ was obtained.
They were sorted, the interval from the 6th to 195th value taken
as representing a~95\%  ``precision interval" characterizing each
method. The results are in the following Table~\ref{table1}, its left part.
It contains also the mean of 200 estimates.

\begin{table}[h]
 %\vspace{5mm}
  \centering
\resizebox{\textwidth}{!}{%
\begin{tabular}{|l|c|c||c|c|}
\hline
 Method & $N=100$: 95\% ``PI'' & Mean & $N=500$: 95\% ``PI'' & Mean \\
 \hline
 Miller   & $(-0.9740,  -0.8410)$ & $-0.9164$ & $(-0.9990,  -0.9550)$ &$ -0.9774$  \\
 BJ       & $(-1.0188,  -0.9189)$ & $-0.9735$ & $(-1.0220,  -0.9739) $& $-0.9987$ \\
 KSVR  & $(-2.2786,    0.0590)$ & $-0.6009$ & $(-1.2042,  -0.5121)$& $-0.7991$\\
 MLE    & $(-1.0400,  -0.9650)$ & $-1.0002 $& $( -1.0100,  -0.9800)$&  $-0.9981$\\
 AFT    & $(-1.0434,  -0.9510)$ & $-0.9989$ & $( -1.0206,  -0.9819)$&  $-1.0016$\\
 \hline
 \end{tabular}}
\caption{Results of estimation of regression parameter $\beta$,
200 repetitions, from data of size $N=100$ and $N=500$.}%
\label{table1}%
 \end{table}

The right part of Table~\ref{table1} shows the results of estimation from
larger data-sets, with $N=500$, each experiment was again repeated
200 times. It is seen that the parametric MLE method and the
estimation in the framework of semi-parametric AFT model have
provided the most reliable estimates of the regression parameter.
The performance of the Buckley and James method was just slightly
weaker, while the method of Miller seemed to be biased
systematically. Simultaneously, all methods have shown a~tendency
to improve its precision with growing data size. The method of
Koul et al., however, has again demonstrated that the large
variability of ``synthetic'' data can cause an instability of
trend estimation.

\section{Concluding remarks}

We have recalled three ``historical'' methods analyzing the LRM
with censored data, compared them mutually, with a~direct
(numerical) solution via the MLE under assumption of normally
distributed errors, and also with a~solution in the framework of
the semi-parametric AFT regression model. The results seem to
support a~conjecture that the BJ estimator is the most reliable
from the three traditional approaches. Very good results were
achieved by the last two methods. It could be said that they
represent contemporary approaches to the problem of linear
regression with censored data. A~variant methods are nowadays
based on the Bayes approach and the use of the Markov chain Monte
Carlo (MCMC) procedures (as for example in Volf and Timkov\'a,
2014).
\cite{VoTi}%

Let us here also mention a~connection of the LR problem with the
Cox model, at least with its special case. Let us assume that the
baseline distribution is the Weibull one, with the distribution
function $F_0(t)=1-\exp\{-(t/a)^b\}$ defined for $t>0$ and with
positive parameters, scale $a$ and shape $b$. Further, let the
dependence on a~covariate $x$ follows the Cox (proportional
hazard) regression model, let us denote the Cox regression
parameter by $\gamma$. Then the cumulative hazard rate equals:
  $$
H(t|X=x)=H_0(t)\cdot \exp(\gamma\,x)= (t/a)^b\cdot
\exp(\gamma\,x)=\left({t\over a(x)}\right)^b, \eqno{(12)}
 $$
where $a(x)=a\cdot \exp\{-\gamma\,x/b\}$. Thus, the scale
parameter $a(x)$ depends on covariate, while the shape parameter
remains $b$. Simultaneously, $H(t|X=x)$ can be rewritten as
 $$
H(t|X=x)=\left({t\cdot \exp(\gamma/b\cdot x)\over a}\right)^b=H_0(t\cdot \exp
(\beta\,x)), \eqno{(13)}
 $$
which corresponds to the AFT model with regression parameter
$\beta=\gamma/b$. Hence, if in the case of Weibull baseline
distribution the estimation uses the Cox model framework, we then
obtain the estimate of $\gamma=\beta\cdot b$ instead of~$\beta$,
while the estimates of $H_0(t)$ in both Cox and AFT models should
be close one to each other (because both are consistent, both
should tend to the Weibull$(a,b)$ cumulative hazard rate).

%%%%%%%%%%%%%
%\newpage

%\renewcommand{\refname}{Literatura}
\begin{thebibliography}{99}
\setlength\itemsep{1pt}

\bibitem{BaNi} Bagdonavi\v{c}ius, V., Nikulin, M. (2002): {\it Accelerated Life
Models, Modeling and Statistical Analysis.}
Chapman \& Hall/CRC, Boca Raton.

\bibitem{BuJa} Buckley, J., James, I. (1979):
Linear regression with censored data. {\it Biometrika} {\bf 66},
429--436.

\bibitem{Cox} Cox, D.\,R. (1972):
Regression models and life-tables (with Discussion). {\it J. R. Stat. Soc.}
{\bf B 34}, 187--220.

\bibitem{JaSm} James, I.\,R., Smith, P.\,J. (1984):
Consistency results for linear regression with censored data. {\it
Ann. Statist.} {\bf 12}, 590--600.

\bibitem{Jin} Jin, Z., Lin, D.\,Y., Ying, Z. (2006):
On least-squares regression with censored data. {\it Biometrika}
{\bf 93}, 147--161.

\bibitem{KaPr}  Kalbfleisch, J.\,D., Prentice, R.\,L. (2002): {\it The
Statistical Analysis of Failure Time Data.} Wiley, New York.

\bibitem{KSVR} Koul, H., Susarla, V., Van Ryzin, J. (1981):
Regression analysis with randomly right-censored data. {\it Ann.
Statist.} {\bf 9}, 1276--88.

\bibitem{LiWe} Lin, D.\,Y., Wei, L.\,J., Ying, Z. (1998): Accelerated
failure time models for counting processes.
{\it Biometrika} {\bf 85}, 605--618.

\bibitem{Mill} Miller, R.\,G. (1976): Least squares regression
with censored data. {\it Biometrika} {\bf 63}, 449--464.

\bibitem{Nov} Nov\'ak, P. (2013): Goodness-of-fit test for the
accelerated failure time model based on martingale residuals. {\it
Kybernetika} {\bf 49}, 40--59.

\bibitem{VoTi} Volf, P., Timkov\'{a}, J. (2014): On selection of optimal
stochastic model for accelerated life testing. {\it Reliability
Engineering \& System Safety} {\bf 131}, 291--297.

\end{thebibliography}

