lifemodels¶

What is it?¶

lifemodels is a Python package for fitting lifespan data to a number of models using MLE. It also implements GLM (generalized linear model) for those lifespan distributions.

It is the next version of deday. deday stands for Demography Data Analyses. The older versions can be found here.

Read the minimal example for a simple workflow in IPython Notebook

Where to find support?¶

Read the documentation for instructions.

If you use deday for published research, please contact me as lifemodels has not yet been published. To install, run pip install git+https://github.com/ctzhu/lifemodels

• Time-to-event data¶

• Engineering: Reliability (MTTF)

• Marketing: Attribution (web cookie expiration, time-dependent binary responce)

• Biomedical science: lifespan, survivorship.

• Censor and missing data¶

• Event time is $[T0, T1)$

• T0==T1, or only know T0, or only know T1, ...

• Time-dependent and Time-independent components:¶

• Benjamin Gompertz & William Makeham.

• Gompertz Makeham law of mortality: $H(t)=\alpha e^{\beta t} + C$

• Example: Estimated probability of a person dying at each age, for the U.S. in 2003

Basic usage¶

· Read the data and create a survival object using .survt_df() method¶

· Data taken from Clive McCay et.al :¶

The effect of retarded growth upon the length of life span and upon the ultimate body size. J Nutr. 1935;10:63–79

· Credited with the discovery of dietray restriction (Eat less live longer)¶

· Two treatments: Diet and Sex. III type diet being the richest.¶

· Death events of rat recorded over ~1500 days.¶

· Specify the model and create the model fit object using .distfit_df() method¶

· Hazard functions of 5 models in this example are:¶

$H(t)=C$

Exponential Decay

$H(t)=(\kappa/\lambda)\cdot(t/\lambda)^{\kappa-1}$

$log(H(t))=log(\kappa \cdot \lambda^{-\kappa})+(\kappa-1)\cdot log(x)$

Weibull model (Mortality increase with time in a power function relationship) 'wb2'

$H(t)=\alpha e^{\beta t}$

$log(H(t))=log(\alpha)+ \beta t$

Gompertz model (Mortality increase with time Exponentially) 'gp2'

$H(t)=\alpha e^{\beta t} - \alpha + C$

Gompertz-Makehamm model 'gp3'

$H(t)=(\kappa/\lambda)\cdot(t/\lambda)^{\kappa-1} + C$

Weibull-Makehamm model 'wb3'

$H(t)=\frac{\alpha e^{\beta t}}{d e^{\beta t} + 1} - \frac{\alpha}{d+1} + C$

Gompertz-logistic model 'gl4'

· MLE parameter estimates¶

· Variance of the parameter estimates¶

It might not always be there, due to float overflow and underflow problems

· log quantile function (time-to-reach-%-of-death estimate): .logq() method¶

.logq(key, q), key is genotype-food combination here. q is that target % of death, default to 50%, therefore by default it returns the estimate for medium lifespan.

.logq() method and .logq_nc() method are similar: the latter ignores the Makehamm term ($C$), if there is one.

· Get the estimated log mean lifespan using .logm() method¶

· ... and the $\sigma^2$ for calculating the confidence interval for survival time using .logt_var()¶

· Plot harzard plot, survival plot and other plots using .plot() method¶

· Plot multidimensional confidence regions with .contourf_2D() method¶

· AKA. Why Gompertz distribution is rarely worked on.¶

· Weibull model is a well-conditionned quadratic function¶

· Gompertz model is an ill-conditionned very non-quadratic function.¶

· Challenge of Gompertz model: precision lost due to overflow/underflow.¶

· Gompertz-Makeham model: multidimensional optimization problem is even harder¶

Experimental Functionality¶

· GLM (generalized linear models) for scale, shape, location, $\alpha$, $\beta$, $\kappa$, $\lambda$... (you name it)¶

· ... and partially fixed models too.¶

· With Weibull $\lambda$ fixed, it is a form of parametric proportional hazard model¶

· With Weibull $\kappa$ fixed, it is a from of ALT (Accelerated life testing) or parametric survival regression model¶

· R does it too, let's compare the solutions¶

R results:¶

> S1 <- survreg(Surv(Time0) ~ Diet*Sex, dist='weibull', data=df)
> summary(S1)

Call:
survreg(formula = Surv(Time0) ~ Diet * Sex, data = df, dist = "weibull")
                  Value Std. Error       z        p
(Intercept)      6.7515     0.0827  81.680 0.00e+00
DietII           0.0776     0.1156   0.671 5.02e-01
DietIII          0.0462     0.1211   0.382 7.03e-01
SexMale         -0.4243     0.1323  -3.208 1.34e-03
DietII:SexMale   0.4459     0.1885   2.366 1.80e-02
DietIII:SexMale  0.4725     0.1881   2.512 1.20e-02
Log(scale)      -0.9506     0.0834 -11.393 4.56e-30

Scale= 0.387 

Weibull distribution
Loglik(model)= -764.8   Loglik(intercept only)= -771.7
    Chisq= 13.89 on 5 degrees of freedom, p= 0.016 
Number of Newton-Raphson Iterations: 7 
n= 106 

· Model both $\lambda$ and $\kappa$¶

· It is an inference problem, biologists like to assign parameters with biological meanings and like to know if any treatment specifically affects one type of biological process but not the others.¶

Notes on Implementation:¶

· log-likelihood calculation is costy¶

· ... and I want to handle interval censored, right censored, mix censored, etc in one shot¶

· ... and in some cases I want to avoid precision lost¶

· Let's do it in FORTRAN¶

mortdist.f90

subroutine lg3LLK(sx1, sx2, sy, su, ss, sc, n)
            double precision, intent(in) :: sx1(n), sx2(n)
            double precision, intent(out) :: sy
            double precision :: su, ss, sc
            integer n, i

            real(kind=10) :: x1(n), x2(n)
            real(kind=10) :: y
            real(kind=10) :: v1, v2, u, s, c

            x1= sx1
            x2= sx2
            u = su
            s = ss
            c = sc
            y = 0.0_10

            do 100 i=1, n
                v1 = exp((u-x1(i))/s)
                if (x1(i) == x2(i)) then
                    y = y+log(v1)-log(1.0_10+v1) &
                         -c*x1(i) &
                         +log(1.0_10/s/(1.0_10+v1)+c)
                else
                    v2 = exp((u-x2(i))/s)
                    y = y+log(v1/(1.0_10+v1)/exp(c*x1(i)) &
                             -v2/(1.0_10+v2)/exp(c*x2(i)))
                end if
    100     continue
            sy = dble(y)
            end subroutine lg3LLK

· Grid-search (to get initial guess) is costy, calculating derivative is costy¶

· scipy.stats build-in continuous distributions are slow¶

· Do them in FORTRAN as well¶

Getting CDF, PDF, SF, and H from one call to avoid calculate the same thing more than once.

mortdist_cpsh.f90

subroutine gp2_cpsh(p, x, cdf, pdf, sf, hzd, m, n)
            double precision, intent(in) :: p(m), x(n)
            double precision, intent(out):: cdf(n), pdf(n), sf(n), hzd(n)
            double precision :: a, b
            double precision :: ebx(n)
            integer :: m, n
            a = p(1)
            b = p(2)
            ebx = exp(b*x)
            hzd = a*ebx
            sf  = exp(-a/b*(ebx-1.0d0))
            pdf = hzd*sf
            cdf = 1.0d0-sf
            end subroutine gp2_cpsh