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
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Reply via email to