On 15/11/12 21:22, David Winsemius wrote:
On Nov 15, 2012, at 5:38 AM, Matthias Ziehm wrote:

Hi all,

Sorry if this has been answered already, but I couldn't find it in the archives 
or general internet.

A Markmail/Rhelp search on:  gompertz survreg  ...brings this link to a reply 
by Terry Therneau. Seems to address everything you asked and a bit more

http://markmail.org/search/?q=list%3Aorg.r-project.r-help+gompertz+survreg#query:list%3Aorg.r-project.r-help%20gompertz%20survreg+page:1+mid:6xdlsmo272oa7zkw+state:results

(Depending on how your mailer breaks URLs you may need to paste it back 
together.)

Thanks! I hadn't read this thread because of the missleading title. However, now I've got follow up questions, on that explanation:

On Feb 5, 2010 at 8:48:08 am, Terry Therneau wrote:
Subject: Re: [R] Using coxph with Gompertz-distributed survival data.
...
the note below talks about how to do so approximately with survreg.
It's a note to myself of something to add to the survival package
documentation, not yet done, and to my embarassment the file has a
time stamp in 1996. Ah well.

Terry Therneau

My document "A Package for Survival Analysis in S" contains statements
about how to fit Gompertz and Rayleigh distributions with the survreg
routine. Nicholas Brouard, in a recent query to this group, quite
correctly states that "Therneau's documentation is a little elliptic for
people not so familiar with extreme value theory".

I've spent the last day trying to work out concrete examples of the
fits. Let me start by saying that I now think my paper's remarks were
overly optimistic. This note will try to indicate why. I will use some
"TeX" notation below: \alpha, \beta, etc for Greek letters.

Hazard functions:

Weibull: p*(\lambda)^p * t^(p-1)
Extreme value: (1/ \sigma) * exp( (t- \eta)/ \sigma)
Rayleigh: a + bt
Gompertz: b * c^t
Makeham: a + b* c^t

The Makeham hazard seems to fit human mortality experience beyond
infancy quite well, where "a" is a constant mortality which is
independent of the health of the subject (accidents, homicide, etc) and
the second term models the Gompertz assumption that "the average
exhaustion of a man's power to avoid death to is such that at theendof
equal infinitely small itervals of time he lost equal portions of his
remaining power to oppose destruction which he had at the commencement
of these intervals". For older ages "a" is a neglible portion of the
death rate and the Gompertz model holds.

The fitting routine depends on the decomposition Y = \eta + \sigma W,
where \eta = \beta_0 + \beta_1 * X_1 + \beta_2 * X_2 + ... is the fitted
linear predictor and W is a distribution in "standard" form. For
instance, if the response time t is Weibull, then y = log(t) follows
this with \eta = log(\lambda) \sigma = 1/p

Clearly

1. The Wiebull distribution with p=2 (sigma=.5) is the same as a
Rayleigh distribution with a=0. It is not, however, the most general
form of a Rayleigh.

2. The (least) extreme value and Gompertz distributions have the same
hazard function, with \sigma = 1/ log(c), and exp(-\eta/ \sigma) = b.
If I do the math correctly from the above given extreme value hazard to Gompertz hazard. It needs to b = 1/ \sigma * exp(-\eta/ \sigma)
Other wise the 1/ \sigma of the Extreme value hazard is missing, isn't it?


It would appear that the Gompertz can be fit with an identity link
function combined with the extreme value distribution. However, this
ignores a boundary restriction. If f(x; \eta, \sigma) is the extreme
value distribution with paramters \eta and \sigma, then the definition
of the Gompertz densitiy is g(x; \eta, \sigma) = 0 x< 0 g(x; \eta,
\sigma) = c f(x; \eta, \sigma) x>=0

where c= exp(exp(-\eta / \sigma)) is the necessary constant so that g
integrates to 1.
Here I got really lost were the addition double exp suddenly comes from and how it fits in.
Given the above I would have thought that:
g(x,b,c) = f(x, \eta=-1/log(c)*log(b*1/log(c)), \sigma= 1/log(c) ) for x>=0

Can anyone clarify these thinga to me, please?

Matthias

If \eta / \sigma is far from 1, then the correction
term will be minimal and survreg should give a reasonable answer. If
not, the distribution can't be fit, nor can it be made to easily conform
to the general fitting scheme of the program.

The Makeham distribution falls into the gamma family (equation 2.3 of
Kalbfleisch and Prentice, Survival Analysis), but with the same range
restriction problem.

In summary, the Gompertz is a truncated form of the extreme value
distribution (Johnson, Kotz and Blakrishnan, Contiuous Univariate
Distri- butions, section 22.8). If one ignores the truncation, i.e.,
assume that negative time values are possible, then it can be fit with
survreg. My original note seems to have been compounded of 3 errors: the
-1 arises from confusing the maximal extreme distribution (most common
in theory books) with the minimal extreme distribution (used in
survreg), the log() term was a typing mistake, and I never noticed the
range restriction.

This is one of the few topics in the report without a worked example as
part of my test library (the Examples section of the package). The
replacement document, currently in early draft, is intended to have a
worked example for every claim and the code for that example in the
appendix. This will, hopefully, cure any other mistakes of this sort.


Is it possible to implement the gompertz distribution as survreg.distribution 
to use with survreg of the survival library?
I haven't found anything and recent attempts from my side weren't succefull so 
far.

I know that other packages like 'eha' and 'flexsurv' offer functions similar to 
survreg with gompertz support. However, due to the run-time environment were 
this needs to be running in the end, I can't use these packages :(

Same questions for the gompertz-makeham distribution.

Many thanks!

Matthias

______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Reply via email to