Many thanks to Prof. Lumley. You are right, the likelihood function could be
written with builtin gamma distribution functions. Now the MLE works
completely fine. :-)
Best Regards,
Kuan-Ta Chen
----- Original Message -----
From: "Thomas Lumley" <[EMAIL PROTECTED]>
To: "Kuan-Ta Chen" <[EMAIL PROTECTED]>
Cc: <[EMAIL PROTECTED]>
Sent: Tuesday, October 19, 2004 1:37 AM
Subject: Re: [R] Survreg with gamma distribution
> On Tue, 19 Oct 2004, Kuan-Ta Chen wrote:
>
> >
> > One million thanks to Prof. Ripley and Prof. Lumley. I think I now have
more
> > understanding regarding survreg with gamma distribution. But one of my
> > problems is still there: in the text of Lee, Wang (2003), there are two
> > "kinds" of parametric fitting: 1) fitting of survival distributions
(like
> > regular probabillity distribution fitting) 2) regression model fitting
> > (mostly assume an accelerated failure time model). Survreg {survival}
> > provides model fitting of (2). But I still have one problem regarding
(1):
> > try to estimate the parameters of gamma distributions for some data.
>
> There aren't really two separate kinds: 1 is a special case of 2, so
> survreg() can do 1.
>
> > For regular gamma distr. fitting, we could use fitdistr (mass) or use
> > optim()/mle() with log-likelihood composed by dgamma()/pgamma(). But
because
> > the data contains (randomly) censored observations, the log-likelihood
> > function must be modified to include the effect of duration of censored
> > observations.
>
> Yes. The loglikelihood is
> pgamma(x,shape,scale=scale,lower.tail=FALSE,log.p=TRUE)
> for a censored observation and
> dgamma(x,shape,scale=scale,log=TRUE)
> for an uncensored observation. No integration necessary.
>
> You might want to work with log(shape) and log(scale) instead, to avoid
> the boundaries at 0.
>
> eg if your data were in variables "times" and "status"
> ll <-function(logshape,logscale){
> -sum(ifelse(status,
> pgamma(times,exp(logshape),scale=exp(logscale),
> lower.tail=FALSE,log.p=TRUE),
> dgamma(times,exp(logshape),scale=exp(logscale),log=TRUE)
> ))
> }
>
> This works in mle() without too much sensitivity to starting values.
>
> -thomas
>
______________________________________________
[EMAIL PROTECTED] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
