[R] Survreg with gamma distribution

[Kuan-Ta Chen]

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" <tlumley at u.washington.edu>
To: "Kuan-Ta Chen" <kuan at ilife.cx>
Cc: <r-help at stat.math.ethz.ch>
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
>