[R] Survreg with gamma distribution
Kuan-Ta Chen
kuan at ilife.cx
Mon Oct 18 19:01:41 CEST 2004
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.
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. To clarify, I've excerpted the log-likelihood function and two
equations of gamma and lambda by taking the first derivation. But
unfortunately, but the loglik function and equations contain integrations
and I can't analytically eliminate them. That's the reason why I used
integration in optim() and always got errors (since I don't have clues to
handle divergent integration.)
The excerpt is from Lee, Wang (2003) p.193 (sorry I don't have another way
to show the complicated equations):
http://kuan.ilife.cx/gammamle.jpg The authors suggest using numerical method
to solve the equation and I don't have any idea to eliminate the
integrations from these equations before optim(). Please give me some hint,
thanks.
Lee, Wang (2003): Elisa T. Lee, John Wenyu Wang, Statistical Methods for
Survival Data Analysis, 3rd edition, 2003
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: Monday, October 18, 2004 11:48 PM
Subject: Re: [R] Survreg with gamma distribution
> On Sun, 17 Oct 2004, Kuan-Ta Chen wrote:
>
> > Hi, all:
> >
> > I find survreg {survival} has provided many distributions such as
weibull,
> > lognormal, etc. But I wonder why it doesn't have the support for gamma
> > distribution since it should be a good distr. in lifetime analysis. Can
> > anybody figure out the reason?
>
> Presumably the actual reason is because Terry Therneau didn't need to use
> the Gamma model. However, all the distributions in survreg are
> location-scale families, which the Gamma is not, so the basic algorithm
> would have to be different.
>
> > I've tried to implement the likelihood function of progressively
censored
> > data for gamma distr. and use optim() to solve the paramemters. The
> > log-likelihood function L contains some integrations.
>
> It shouldn't have to: we do have pgamma() built in (and digamma, trigamma,
> etc for derivatives).
>
> > I use tryCatch() to
> > capture the error when integration lead to divergence and return Inf.
> > But if consequent two calls to the objective function return Inf,
optim()
> > will raise errors:
> >
> > Error in optim(c(ga, 1/la), fr, method = "BFGS") :
> > non-finite finite-difference value [1]
> >
> > What can I do except for choosing better initial values?
>
> Choose better initial values. You should be able to get quite good
> initial values for regression coefficients by using survreg on a lognormal
> distribution, since Gamma and lognormal distributions agree pretty well
> except in the extreme tails. You could then try getting the shape
> parameter by matching the variance of the Gamma to the variance of the
> fitted lognormal.
>
> > The last question, by its name "survreg", survreg does its job by
> > regression,
> > but why p.75 in Tableman, Kim (2004) said that "We use the S function
> > survReg to fit parametric models (with the MLE approach)...". Does
survreg
> > use regression or MLE approach?
>
> Its job *is* regression. It uses maximum likelihood to fit a regression
> model.
>
> -thomas
>
More information about the R-help
mailing list