2
$\begingroup$

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?

$\endgroup$
1
$\begingroup$

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))$$

The mean lifetime would be

$$\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)

#Define your likelihood function:
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]))))))
}

#Your data:
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))

#Your fit:
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)

I have added the survival plot below for your reference. Given your data, the mean lifetime would be:

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

km

$\endgroup$
  • 1
    $\begingroup$ For what it's worth, the OP asked for an EM algorithm, while this is using quasi-Newton's method. With that said, I'm pretty confident that this will be faster than actually using the EM algorithm. $\endgroup$ – Cliff AB Jul 12 '19 at 5:01
  • $\begingroup$ That's a fair point. I glossed over the EM part as it seemed as though the OP thought that the EM algorithm was required to fit the distribution. Obviously it isn't, as we can just directly maximize the likelihood directly using numerical methods. $\endgroup$ – Ed P Jul 12 '19 at 5:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.