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" <[EMAIL PROTECTED]> To: "Kuan-Ta Chen" <[EMAIL PROTECTED]> Cc: <[EMAIL PROTECTED]> 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 > ______________________________________________ [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