[R] Estimation parameters of lognormal censored data
R. Michael Weylandt
michael.weylandt at gmail.com
Wed Aug 29 16:33:37 CEST 2012
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)
>
>
>
>
More information about the R-help
mailing list