Hi, I am trying to get the maximum likelihood estimator for lognormal distribution with censored data;when we have left, interval and right censord. I built my code in R, by writing the deriving of log likelihood function and using newton raphson method but my estimators were too high " overestimation", where the values exceed the 1000 in some runing of my code. is there any one can tell me where the error can be and help me in how I can get the estimators with right way?. kind regards [[alternative HTML version deleted]]
R. Michael Weylandt
2012-Aug-29 04:21 UTC
[R] Estimation parameters of lognormal censored data
On Tue, Aug 28, 2012 at 9:58 PM, Salma Wafi <salmawafi76 at yahoo.com> wrote:> Hi, I am trying to get the maximum likelihood estimator for lognormal distribution with censored data;when we have left, interval and right censord. I built my code in R, by writing the deriving of log likelihood function and using newton raphson method but my estimators were too high " overestimation", where the values exceed the 1000 in some runing of my code. > is there any one can tell me where the error can be and help me in how I can get the estimators with right way?. > kind regardsThis sounds like a rather complicated calculation and would be much easier to follow along with if you could provide a reproducible example: this site provides many useful pointers. http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example Cheers, Michael> [[alternative HTML version deleted]] > > > ______________________________________________ > R-help at 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. >
R. Michael Weylandt
2012-Aug-29 14:33 UTC
[R] Estimation parameters of lognormal censored data
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) > > > >