Hazard function of a gamma distribution

Question

The system we are working on is biological, more specifically the distribution of specific events across a chromosome. This can be thought of as 1D array (the chromosome) across which points can be chosen (event positions). We have mapped the positions of these events experimentally and then calculated the inter-point distances (IPDs) (the distances between each successive event) - see below.

enter image description here

We have fitted a gamma-distribution to the calculated IPDs:

enter image description here

Gamma distribution

a = 2.7863
b = 39498.2

Which can be reconstructed separately using:

dist = gampdf(1:450000,2.7863,39498.2);

Question 1: I now want to convert this gamma distribution to a hazard function. Am I right in thinking that this will essentially plot out the probability of an event occurring at position x relative to an event at position 0?

To expand upon that, we have previously observed that our IPDs do not match those expected if event distributions across chromosomes were totally random. We are now working under the hypothesis that each event that occurs, creates a zone around it which makes it less likely another will occur nearby. By using a hazard function, we want to obtain measures of this reduction in probability at different distances.

Question 2: Could anybody provide the hazard function from this distribution so I know what I am aiming for? This is because I am having difficulty finding a way to do this on my chosen software and need to be able to check my output.

IPD data used to produce the gamma-distribution:

user78229 · Accepted Answer · 2015-05-27 13:47:08Z

2

First point to note is that your a-value of 2.7863 suggests that the rate of these events is increasing over time. Given that, your assumption of randomness may not hold.

Paul Allison's book Survival Analysis Using SAS, on page 295 has a formula for calculating the hazard for the gamma distribution that is independent of the software. It's a bit hairy:

k=1/(shape parameter*shape parameter)
do t=min to max by inc  (min, max are obs'd values of time, incremented)
  u=(t*exp(-m)**(1/scale parameter)
  f=abs(shape)*(k*u**shape)**k*exp(-k*u**shape)/(s*gamma(k)*t)
  Surv=1-probofgamma(k*u**shape,k)
    if shape<0 then Surv=1-Surv
  Hazard=f/Surv
  output
end

where m = mean of the predicted values from your gamma dist survival model

Hope this helps...

answered May 27, 2015 at 13:47

user78229

10.9k2 gold badges25 silver badges46 bronze badges

3

$\begingroup$ 1. Why does this particular value of the Gamma shape parameter $a$ indicate something is "increasing over time"? What exactly is your criterion? 2. That looks like SAS code: it's hardly "independent of the software." It would be preferable to provide the mathematical formula. $\endgroup$
– whuber ♦
Commented May 27, 2015 at 18:05

Add a comment |

Ytsen de Boer · Accepted Answer · 2016-09-20 22:59:12Z

Question 1

No, not really. The hazard function plots out a number, proportional to the probability that you find the next event in the interval $\texttt{lim}_{\Delta\,x\rightarrow\infty}[x, x+\Delta\,x]$ as you increase $x$ and have not yet come across the next event.

In general I don't believe it is a probability distribution.

Question 2

The hazard function $h(x)$ for a distribution is defined as the ratio between its probability density function and its survival function.

Given your fit (which looks very good) it seems fair to assume the gamma function indeed.

Both the pdf and survival function can be found on the Wikipedia page of the gamma distribution.

So (check this) I got:

$h(x) = \frac{x^{a-1}\,e^{-x/b}}{b^{a}\left(\Gamma(a) - \gamma(a, x/b)\right)}$

Here $\gamma$ is the lower incomplete gamma function

Beyond your question

If your measurements allow the parameter b to be integer, then the good gamma fit (with positive parameters) also indicates a good Erlang fit. In that case the hazard function may become easier to compute yourself:

$h(x)=\frac{\left(x/a\right)^{b-1}}{a(b-1)!\sum_{i=0}^{b-1}\frac{(x/a)^{i}}{i!}}$

Also, in case of using the Erlang distribution, your IPDs could be thought of as being "built up" out of about $b=40K$ smaller distances, which are, in turn, exponentially distributed, with parameter $a=2.8$.

This may help explain the size of the "zone" you mention.

Tim · Accepted Answer · 2016-09-20 19:52:43Z

1

Hazard function can be defined as

$$ h(x)=\frac{f(x)}{1-F(x)}=\frac{f(x)}{S(x)} $$

where $f$ is probability density function and $F$ is cumulative distribution function ($S$ is survival function). So to calculate hazard function you basically need pdf and cdf. This is what code provided by @MikeHunter does. Check here for pdf, cdf and hazard function formulas for gamma distribution.

answered Sep 20, 2016 at 19:52

Tim

141k26 gold badges270 silver badges512 bronze badges

Add a comment |

Stack Exchange Network

Hazard function of a gamma distribution

3 Answers 3

Question 1

Question 2

Beyond your question

Your Answer

Not the answer you're looking for? Browse other questions tagged
distributions
matlab
survival
gamma-distribution
or ask your own question.

Hot Network Questions

Hazard function of a gamma distribution

3 Answers 3

Question 1

Question 2

Beyond your question

Your Answer

Sign up or log in

Post as a guest

Not the answer you're looking for? Browse other questions tagged distributionsmatlabsurvivalgamma-distribution or ask your own question.

Related

Hot Network Questions

Not the answer you're looking for? Browse other questions tagged
distributions
matlab
survival
gamma-distribution
or ask your own question.