Inverse Gaussian Frailty

Why the Frailty Distribution Matters

A frailty model describes individual event risk as a conditional hazard multiplied by a positive, unobserved risk factor. The distribution chosen for that factor determines how conditional survival averages across individuals. For inverse Gaussian frailty, this averaging can be done with a closed-form Laplace transform, giving explicit expressions for marginal survival and hazard.

The formulas below concern a single event time conditional on fixed covariates. They are also building blocks for recurrent-event and clustered-event models. The key distinction is between a hazard conditional on an individual's frailty and the hazard after that frailty has been averaged over the population.

Definition: A mean-one inverse Gaussian frailty $Z$ with variance parameter $\theta>0$, written $Z\sim\operatorname{IG}(1,1/\theta)$, has density $$ f_Z(z)=\frac{1}{\sqrt{2\pi\theta z^3}}\exp\left\{-\frac{(z-1)^2}{2\theta z}\right\},\qquad z>0. $$ This parameterization gives $E(Z)=1$ and $\operatorname{Var}(Z)=\theta$. As $\theta$ tends to zero, the distribution converges to a point mass at $Z=1$, the no-frailty limit.

Suppose that, conditional on frailty $Z=z$ and covariates $x$, the event hazard at time $t$ is $$ h(t\mid Z=z,x)=z h_0(t)\exp(x^\mathsf{T}\beta), $$ where $h_0(t)$ is a baseline hazard and $\beta$ is a vector of regression coefficients. Define the cumulative baseline hazard $H_0(t)=\int_0^t h_0(u)\,du$ and $$ A(t,x)=H_0(t)\exp(x^\mathsf{T}\beta). $$ Assume $h_0(t)\geq 0$, the cumulative hazard exists over the time interval under consideration, $\theta>0$, and the frailty is independent of the covariates for the population being described. Conditional on $Z$ and $x$, the survival function is $\exp\{-ZA(t,x)\}$. No censoring assumption is needed for this mathematical identity; assumptions about censoring matter when estimating the model from data.

Deriving the Laplace Transform

The Laplace transform of a nonnegative random variable $Z$ is $L_Z(s)=E\{\exp(-sZ)\}$, for $s\geq0$. For inverse Gaussian frailty, substitute its density into this expectation. Expanding the quadratic term in the density gives $$ L_Z(s)=\frac{\exp(1/\theta)}{\sqrt{2\pi\theta}}\int_0^\infty z^{-3/2}\exp\left\{-\left(s+\frac{1}{2\theta}\right)z-\frac{1}{2\theta z}\right\}\,dz. $$

For positive $a,b$, the integral used here is $$ \int_0^\infty z^{-3/2}\exp(-az-b/z)\,dz=\sqrt{\frac{\pi}{b}}\exp(-2\sqrt{ab}). $$ To justify it, set $u=\sqrt{b/z}$. The integral becomes $2/\sqrt b$ times $J(\sqrt{ab})$, where $J(c)=\int_0^\infty \exp(-u^2-c^2/u^2)\,du$. For $c>0$, differentiation under the integral and the substitution $v=c/u$ give $J'(c)=-2J(c)$. Since $J(0)=\sqrt{\pi}/2$, solving this differential equation gives $J(c)=(\sqrt{\pi}/2)\exp(-2c)$, establishing the identity.

In the transform integral, take $a=s+1/(2\theta)$ and $b=1/(2\theta)$. The factor $\sqrt{\pi/b}$ cancels the density's leading constant, apart from $\exp(1/\theta)$. Simplifying the exponent yields the transform:

Result: The mean-one inverse Gaussian frailty has Laplace transform $$ L_Z(s)=\exp\left\{\frac{1-\sqrt{1+2\theta s}}{\theta}\right\},\qquad s\geq0. $$ At $s=0$, $L_Z(0)=1$, as required for a probability transform.

The transform also checks the mean and variance implied by the parameterization. Differentiating gives $$ L_Z'(s)=-\frac{L_Z(s)}{\sqrt{1+2\theta s}},\qquad L_Z''(0)=1+\theta. $$ For a nonnegative random variable, $L_Z'(0)=-E(Z)$ and $L_Z''(0)=E(Z^2)$. Thus $E(Z)=1$, and $\operatorname{Var}(Z)=E(Z^2)-\{E(Z)\}^2=(1+\theta)-1=\theta$. These derivative checks are useful for detecting a wrong scale parameter or a missing factor of two.

Worked Example: Checking the Mean and Variance

Suppose $\theta=0.25$. Use the transform to verify the first two moments and interpret the variance parameter.

At $s=0$, $L_Z(0)=\exp\{(1-\sqrt{1})/0.25\}=1$. The first derivative at zero is $L_Z'(0)=-1$, so $E(Z)=1$. The second derivative is $L_Z''(0)=1+\theta=1.25$, so $E(Z^2)=1.25$.

Therefore $\operatorname{Var}(Z)=1.25-1^2=0.25$, and the standard deviation is $\sqrt{0.25}=0.5$. Frailty is dimensionless, so these moments have no time units. A value of $\theta=0.25$ describes variation in individual risk around a population mean multiplier of one; it does not mean that every person's multiplier equals one.

Marginal Survival and Hazard

Marginal survival averages conditional survival over the frailty distribution. In symbols, for a fixed covariate value $x$, $$ S_m(t\mid x)=E_Z\{S(t\mid Z,x)\}=E_Z\{\exp[-ZA(t,x)]\}=L_Z\{A(t,x)\}. $$ Substituting the inverse Gaussian transform gives

S_m(t\mid x)=\exp\left\{\frac{1-\sqrt{1+2\theta A(t,x)}}{\theta}\right\}.

The marginal hazard is defined as the negative time derivative of the log marginal survival. Since $A'(t,x)=h_0(t)\exp(x^\mathsf{T}\beta)$, differentiation gives

h_m(t\mid x)=-\frac{d}{dt}\log S_m(t\mid x) =\frac{A'(t,x)}{\sqrt{1+2\theta A(t,x)}} =\frac{h_0(t)\exp(x^\mathsf{T}\beta)} {\sqrt{1+2\theta H_0(t)\exp(x^\mathsf{T}\beta)}}.

The denominator is dimensionless, so the marginal hazard has the same units as $h_0(t)$, usually events per unit time. The marginal hazard is not obtained by simply replacing $Z$ with its mean in the conditional hazard: that would give $h_0(t)\exp(x^\mathsf{T}\beta)$ at every time and would omit the denominator. Although $E(Z)=1$, survivors at later times are a selected subset of the original population.

Worked Example: Computing Marginal Survival

Consider a hypothetical individual with $\theta=0.4$ and $A(t,x)=1.2$ at a specified time. Calculate the marginal survival probability at that time.

Using the transform, $$ S_m(t\mid x)=\exp\left\{\frac{1-\sqrt{1+2(0.4)(1.2)}}{0.4}\right\}. $$ The quantity inside the square root is $1+0.96=1.96$, so its square root is $1.4$. The exponent is $(1-1.4)/0.4=-1$, and therefore $S_m(t\mid x)=\exp(-1)\approx0.3679$.

This is the model-implied population survival at that time, after averaging over frailty. For comparison, survival at the mean frailty $Z=1$ is $\exp(-1.2)\approx0.3010$. These are different quantities: the second is conditional survival for an individual with frailty exactly one, not the population average.

Worked Example: Computing the Marginal Hazard

Use the same hypothetical setting as above, and suppose $h_0(t)=0.30$ per year and $\exp(x^\mathsf{T}\beta)=1.5$ at the time of interest. The conditional hazard for $Z=1$ is $0.30(1.5)=0.45$ per year. Find the marginal hazard.

Here $A'(t,x)=0.45$ per year and $A(t,x)=1.2$. Thus $$ h_m(t\mid x)=\frac{0.45}{\sqrt{1+2(0.4)(1.2)}}=\frac{0.45}{1.4}=0.3214\text{ per year}. $$

The population's marginal hazard is $0.3214$ per year at this time, whereas the conditional hazard at $Z=1$ is $0.45$ per year. The difference follows from averaging over frailties among those still event-free; it is not a change in the units or a new baseline hazard.

Checks, Comparisons, and Limitations

Several checks help establish whether these expressions have been applied coherently. At $t=0$, if $H_0(0)=0$, then $A(0,x)=0$, $S_m(0\mid x)=1$, and the marginal hazard equals the conditional hazard multiplier $h_0(0)\exp(x^\mathsf{T}\beta)$, when defined. For $\theta$ approaching zero, $\sqrt{1+2\theta A}$ approaches one and the transform approaches $\exp(-A)$; the no-frailty survival and hazard are recovered. For $A\geq0$, the transform gives survival between zero and one, and the denominator in the hazard is at least one.

A useful comparison is with gamma frailty, whose mean-one Laplace transform is $(1+\theta s)^{-1/\theta}$ under the same variance convention. The lognormal frailty discussed in the earlier material on Lognormal Frailty and Gaussian Random Effects generally has no elementary closed-form Laplace transform. The inverse Gaussian transform is explicit, but choosing it rather than gamma or lognormal is a modeling decision, not a consequence of algebra. Different frailty distributions can imply different marginal survival and hazard shapes even when their means and variances agree.

These formulas assume the specified inverse Gaussian distribution and the stated multiplicative conditional hazard. If frailty is associated with covariates, a single common marginal expression based on the same frailty distribution may not describe all covariate groups. Similarly, a transform calculation alone does not establish that this distribution fits observed data. When fitting a model, inspect whether estimates are numerically stable, whether $\theta$ is nonnegative and plausible, and whether predicted survival and hazards behave sensibly over the observed time range.

Wrong: Treat $E(Z)=1$ as a reason to replace $Z$ by one before calculating survival. Correction: In general, $E\{\exp(-ZA)\}\neq\exp\{-E(Z)A\}$; average conditional survival using the Laplace transform.
Wrong: Use the marginal survival formula but call the conditional hazard the marginal hazard. Correction: Differentiate the log marginal survival; the resulting denominator depends on cumulative hazard.
Wrong: Put time directly into the transform as $s=t$. Correction: Its argument is $A(t,x)$, the cumulative conditional hazard, which is dimensionless.
Wrong: Interpret $\theta$ as the frailty standard deviation. Correction: Under this convention, $\theta$ is the variance and $\sqrt{\theta}$ is the standard deviation.

Key takeaway: For mean-one inverse Gaussian frailty, first calculate the covariate-specific cumulative conditional hazard $A(t,x)$, then evaluate $L_Z\{A(t,x)\}$ for marginal survival. Differentiate its log to obtain the marginal hazard; do not substitute the mean frailty into conditional survival or hazard and assume that averaging is complete.

Check Your Understanding

Answer each question, then check the brief solutions below.

For $\theta=0.5$, what are $E(Z)$ and $\operatorname{Var}(Z)$?
With $\theta=0.5$ and $A=2$, calculate $L_Z(A)$.
If $A'(t,x)=0.60$ per year, $A(t,x)=1$, and $\theta=1$, calculate the marginal hazard.
Why does $E(Z)=1$ not imply that marginal survival is $\exp(-A)$ when $\theta>0$?

Solutions. (1) The parameterization gives $E(Z)=1$ and $\operatorname{Var}(Z)=0.5$. (2) $L_Z(2)=\exp\{(1-\sqrt{1+2(0.5)(2)})/0.5\}=\exp\{2(1-\sqrt3)\}=\exp(-1.4641)\approx0.2312$. (3) $h_m=0.60/\sqrt{1+2(1)(1)}=0.60/\sqrt3\approx0.3464$ per year. (4) Marginal survival is $E\{\exp(-ZA)\}$, which is generally not $\exp\{-E(Z)A\}$; the transform accounts for frailty variation.

What You'll Learn