# How to find the mean lifetime of right-censored data?

Suppose I have a group of patients (the following are their lifetime, * means this individual is censored): 30,67,79*,82*,95,148,170,171,176,193,200,221,243,261,262,263,399,414,446,446*,464,777

The distribution of the lifetimes for these patients is a log-normal distribution. I'm trying to determine the MLE of mu and sigma by using the EM algorithm, and find the mean lifetime of patients in these group. I'm new to EM algorithm so I have no idea what are the steps of EM algorithm and how to find the mean in a right-censored data. Any hints pls?

There are many existing methods of estimating a right-censored log-normal survival model. If we start with the mathematics behind the model, we can build to implementing a solution in a particular software.

Let's assume we have a random variable $$T$$ which represents the survival time of some group. As you have stated, $$T$$ is log-normally distributed:

$$T\sim LN(\mu,\sigma)$$

or alternatively, $$\log T$$ is normally distributed:

$$\log T\sim N(\mu,\sigma^{2})$$

Now, we also know we have a sample of these lifetimes and whether or not these lifetimes are right-censored, $$\{t_{i},c_{i}\}\,\forall i$$. In order to fit our distribution for lifetimes to our data, we need to formulate a likelihood function:

$$L=\prod_{i=1}^{n}(f(t_{i};\mu,\sigma))^{1-c_{i}}(1-F(t_{i};\mu,\sigma))^{c_{i}}$$

where $$f$$ and $$F$$ denote the density function and cumulative distribution function of $$T$$, respectively. Here we are assuming that a right-censored observation implies $$c_{i}=1$$ and non-censored observation implies $$c_{i}=0$$. The interpretation of this function is that each non-censored observation contributes $$f(t_{i})$$ to the likelihood function whereas the censored observations only contribute $$1-F(t_{i})$$ because this represents the probability that the lifetime is greater than $$t_{i}$$, $$\mathbb{P}(T>t_{i})$$.

Essentially, we need to maximize this likelihood function with respect to the parameters $$(\mu,\sigma)$$. We can simplify as follows:

$$\log L=\sum_{i=1}^{n}(1-c_{i})\cdot\log(f(t_{i};\mu,\sigma))+c_{i}\cdot\log(1-F(t_{i};\mu,\sigma))$$

$$\mathbb{E}[T]=e^{\mu+0.5\sigma^{2}}$$

Maximization of the likelihood function can be achieved numerically in any software you choose. In R this would be implemented as follows:

library(survival)

lik=function(pars,c,d) {
return(-sum((1-c)*log(pmin(10^22,pmax(10^-22,dlnorm(d,pars[1],pars[2]))))+c*log(pmin(10^22,pmax(10^-22,1-plnorm(d,pars[1],pars[2]))))))
}

d=c(30,67,79,82,95,148,170,171,176,193,200,221,243,261,262,263,399,414,446,446,464,777)
c=c(0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0)

#Estimate the parameters:
f=stats::optim(c(5,0.7),
fn=lik,
d=d,
c=c,
lower=c(-10,0.0001),
upper=c(Inf,Inf),
method="L-BFGS-B",
control=list(maxit=10^8))

f

#Visualize the fit to data (not an excellent measure of goodness-of-fit by any means):
hist(d,breaks=20,freq=FALSE,main="",xlab="T")
lines(seq(0,max(d),length.out=10^4),dlnorm(seq(0,max(d),length.out=10^4),f$$par[1],f$$par[2]),col="blue",lwd=2)

#Plot survival curve against Kaplan-Meier curve:
s=Surv(d,as.numeric(!c),type="right",origin=0)
plot(s,xlab="T",ylab="S(t)")
lines(seq(0,max(d),length.out=10^4),1-plnorm(seq(0,max(d),length.out=10^4),f$$par[1],f$$par[2]),col="blue",lwd=2)


> exp(f$$par[1]+0.5*f$$par[2]^2)