[R] fitting lognormal censored data

Fri Aug 31 11:43:57 CEST 2012

On 31-08-2012, at 03:54, Salma Wafi wrote:

> Hi ,
> I am trying to get some estimator based on lognormal distribution when we have left,interval, and right censored data. Since, there is now avalible pakage in R can help me in this, I had to write my own code using Newton Raphson method which requires first and second derivative of log likelihood but my problem after runing the code is the estimators were too high. with this email ,I provide simple example for estimation parameters for lognormal when we have only right censored data, and I tried to estimate the parameters using three methods, Surv pakage, optim function, and Newton Raphson calculation. For the Surv pakage and Optim function, I got similar results of estimation values to the true values, but for the Newton Raphson, I got problem with estimation, where it was  too high" overestimation". Is there any one cen help me to know my problem when can be or recommend me to read some good references that can be useful?
> Best regards.
>  Below is the example,  
>  
> cur=curd=cen=cens=array(1,100) 
> Cur=RealCur=realcensoring=realcured=array(1,20)
> ExpCure=Bias=RealCure=array(1,21)
> ##################################
> Z1=c(rbinom(100,1,0.5))
> Z2=c(rbinom(100,1,0.5))
> dat1<-data.frame(time=rlnorm( 100,2,0.8),Censored=rbinom(100,1,0.9),Cured=rbinom(100,1,.3))
> 
> dat2<-dat1[order(dat1[,1]),]                  # order the data #
> for (i in 1:10)
> {
> dat2$Cured[i+90]=0         #Long term survivors/10 individuals#
> dat2$Censored[i+90]=0 
> dat2$time[i+90]=dat2$time[90]
> }
> cens<-c(dat2$Censored)     #censored status #
> curd<-c(dat2$Cured)        #cured status #
> tim<-c(dat2$time)          #lifetimes #
> X1<-c(Z1)                  #X1 Covariate=Age #
> X2<-c(Z2)                  #X2 Covariate=Type of treatment(1=chemo,0=radio) #
> 
> L1<-length(cens)           #number of censored#
> for (j in 1:L1)
> {
> if ((cens[j]==1)&(curd[j]==0)) {(cen[j]=1)&(cur[j]=1)}
> 
> else {(cen[j]=cens[j])&(cur[j]=curd[j])}
> } 
> Now, the following is my data:
> 
> #######    My Data only with uncensored and right censored  ####################
> data=data.frame(Ti=dat1$time,Censored=cen)
> 
> #################### Estimation using Surv pakage ############################
> library(survival)
> survreg(Surv(Ti, Censored)~1, data=data, dist="lognormal")
> ########### Seperating the data for simply using optim function and  Newton Raphson ####
> 
> dat2<-c(split(data[1:2],data$Censored==1))  # Split the data(cen/uncen) # 
> tun<-c((dat2$'TRUE')$Ti)          # Life times data set for uncensored #
> Stat.uncen<-c((dat2$'TRUE')$Censored) # uncensored Status data set #
> tcen<-c((dat2$'FALSE')$Ti)           # Life times data set for censored #
> Stat.cen<-c((dat2$'FALSE')$Censored)  # censoring Status data set #
> ##################### Estimation using optim function  ################
> ml= function(par){mu=par[1]
> s=par[2]
> -sum(dlnorm(tun,mu,s,log=TRUE))-sum(1-plnorm(tcen,mu,s))}
> 
> Max=optim(c(0.5,0.5),ml)
> param=c(Max$par)
> 
> ## My Problem is when I try to estimate parameters using newton raphson  ###
> 
> ########### intial values ##############
> mu=1
> s=1
> ############# Derivative for mu ##########
> F1= function(par){ mu=par[1]
> s=par[2]
> sum((log(tun)-mu)/s^2)+ sum(((1/sqrt(2*pi))*
> exp(-((log(tcen)-mu)^2)/2*s^2))/(s^2*(1-plnorm(tcen,mu,s))))  }
> ############## Derivative for sigma" s "  #######
> 
> F2= function(par){
> mu=par[1]
> s=par[2]
> sum((log(tun)-mu)^2/s^3 -1/s)+ sum((log(tcen)-mu)*((1/sqrt(2*pi))*
> exp(-((log(tcen)-mu)^2)/2*s^2))/(s^2*(1-plnorm(tcen,mu,s))))}
> ############### Total Function  ########################
> F=function(par){
> F=matrix(0,nrow=2)
> F[1]=F1(par)
> F[2]=F2(par)
> F           }
> ################  the Jacobian matrix, a 2 x 2 matrix  ###############
> d11=d12=d21=d22=array()
> J=function(par){
> j=matrix(0,ncol=2,nrow=2)
> # The format of J is 2 by 2#
> d11=0; d12=0;
> d21=0;d22=0
> ######################## second Derivative for mu ##########
> d11 = function(par){
> mu=par[1]
> s=par[2]
> sum(-1/s^2)-sum((1/s^2)*
> (((1/sqrt(2*pi))*exp(-((log(tcen)-mu)^2)/2*s^2))/s*((1-plnorm(tcen,mu,s))))*
> (-(log(tcen)-mu)/s +((1/sqrt(2*pi))*exp(-((log(tcen)-mu)^2)/2*s^2))/
> (s*(1-plnorm(tcen,mu,s)))))}
> 
> ################ Derivative for mu/s   ##########################################
> d12= function(par){
> mu=par[1]
> s=par[2]
> -sum(2*(log(tun)-mu)/s^3)-sum((1/s^2)*(((1/sqrt(2*pi))*
> exp(-((log(tcen)-mu)^2)/2*s^2))/(s*(1-plnorm(tcen,mu,s)))+((log(tcen)-mu)/s)*
> ((1/sqrt(2*pi))*exp(-((log(tcen)-mu)^2)/2*s^2))/(s*(1-plnorm(tcen,mu,s)))*
> (((1/sqrt(2*pi))*
> exp(-((log(tcen)-mu)^2)/2*s^2))/(s*(1-plnorm(tcen,mu,s)))-(log(tcen)-mu)/s)))}
> d21=d12
> ######## second Derivative for s    ###########################
> d22=function(par){
> mu=par[1]
> s=par[2]
> sum((-3*(log(tun)-mu)^2/s^4)+1/s^2)-sum((1/s^2)*((log(tcen)-mu)/s)*(((1/sqrt(2*pi))*
> exp(-((log(tcen)-mu)^2)/2*s^2))/(s*(1-plnorm(tcen,mu,s)))+((1/sqrt(2*pi))*
> exp(-((log(tcen)-mu)^2)/2*s^2))/(s*(1-plnorm(tcen,mu,s)))+((log(tcen)-mu)/s)*(((1/sqrt(2*pi))*
> exp(-((log(tcen)-mu)^2)/2*s^2))/(s*(1-plnorm(tcen,mu,s))))*(((1/sqrt(2*pi))*
> exp(-((log(tcen)-mu)^2)/2*s^2))/(s*(1-plnorm(tcen,mu,s)))-(log(tcen)-mu)/s)))}
> 
> ############  The R codes for Newton's method #####################################
> 
> j[1,1]=d11(par); j[1,2]=d12(par);
> j[2,1]=d21(par); j[2,2]=d22(par)
> j           }
> out1=outs=array()
> par=c(mu,s)        #### initial values #####
> v=matrix(0,ncol=length(par))
> for (i in 1:30)
> {
> v=solve(J(par),F(par))
> par=par-v                             ############ here see the values are to high ###########
> mu[i+1]=par[1]
> s[i+1]=par[2]
> }
> out1=c(mu)                        ############ the value of parameter at each step , no convergence ####################
> outs=c(s)

Your code is a mess and unreadable.

Assuming that the function ml is what you want to minimize, I'm quite sure that the derivatives you are calculating are quite incorrect.
Did you check them?
Assuming that:

- function F gives the gradient of ml
- function J gives the Hessian of ml

we can check these with package numDeriv.

To make the following reproducible I inserted

set.seed(123)

at the start of your code and at the end of your code I inserted

library(numDeriv)
ml.grad <- function(par) grad(ml,par)    
ml.hessian <- function(par) hessian(ml,par)

pstart <- c(0.5,.5)
ml(pstart)
F(pstart)
ml.grad(pstart)
J(pstart)
ml.hessian(pstart)

Running the code in R gives

> ml(pstart)
[1] 570.405
> F(pstart)
             [,1]
[1,] 1.770876e+15
[2,] 7.247892e+15
> ml.grad(pstart)
[1]  -459.6117 -1423.2860
> J(pstart)
              [,1]          [,2]
[1,] -2.634033e+01 -2.567008e+31
[2,] -2.567008e+31 -2.101266e+32
> ml.hessian(pstart)
          [,1]     [,2]
[1,]  326.3977 1826.317
[2,] 1826.3172 9187.053

Functions F and J are most likely not returning correct values (I have no reason to question the results of numDeriv).
So it is not surprising that your attempt at Newton fails miserably.

You can also run this 

(p1.opt <- optim(pstart,ml))
(p2.opt <- optim(pstart,ml, gr=ml.grad, method="BFGS"))
(p3.opt <- optim(pstart,ml, gr=ml.grad, method="CG"))
(p4.opt <- optim(pstart,ml, gr=ml.grad, method="L-BFGS-B"))
ml.grad(p1.opt$par)
ml.grad(p2.opt$par)
ml.grad(p3.opt$par)
ml.grad(p4.opt$par)

which will show you that optim is perfectly capable of minimizing function ml.
I cannot judge whether the result is compatible with the result of survreg. That's up to you.

Finally, it is often not a good idea to use an equation solver to find a solution to gradient of function=0 as a way of locating a minimum.
If you insist on doing that, have a look at packages nleqslv and BB.

regards,

Berend