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)
[[alternative HTML version deleted]]
______________________________________________
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.
