[R] Survreg with gamma distribution

[Thomas Lumley]

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