Salma, I don't know much about survival modelling, but many smart folks do, so I'm forwarding this back to R-help so it gets a wider audience. In fact, we generally like to keep replies "on list" 1) so that they are archived for future googlers; 2) so that you have access to a wider pool of respondents who can help you across time zones and 3) yell and scream when I am wrong (which I often am!) Cheers, Michael
On Wed, Aug 29, 2012 at 3:05 AM, Salma Wafi <salmawafi76 at yahoo.com> wrote:
Dear Michael,
Thanks for your response and help. Dear, actually, I am trying to get
estimator for cure fraction 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". I hope in
this time my English was good to understand. Below is the example, it is
also attached with this email.
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])}
}
####### My Data only with uncensored and right censored
####################
data=data.frame(Ti=dat1$time,Censored=cen)
#################### using Surv pakage ############################
library(survival)
survreg(Surv(Ti, Censored)~1, data=data, dist="lognormal")
########### Seperating the data for simply using optim 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 #
##################### using optim ################
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)