Weibull analysis

The Weibull distribution is extensively used in survival and reliability analysis to model the waiting time for an event to occur. It depends on two positive parameters a and b,

load("distrib")$
pdf_weibull(t, a, b);

$${{a\,\left({{t}\over{b}}\right)^{a-1}\,e^ {- \left({{t}\over{b}} \right)^{a} }}\over{b}}$$

This is the aspect of the density function for different values of a, keeping $b=2$.

load("draw")$
b: 2 $
fpprintprec: 4 $
draw(
  terminal   = animated_gif,
  delay      = 20,
  dimensions = [400, 300],
  makelist(
    gr2d(
      title            = "Weibull density function",
      background_color = orange,
      grid             = true,
      yrange           = [-0.2, 2],
      label_alignment  = left,
      label([string('a = a), 2.5, 1.4]),
      label([string('b = 2.0), 2.5, 1.1]),
      color            = black,
      explicit(pdf_weibull(t, a, b), t, 0, 5)),
    a, 0.2 , 6, 0.2)   ) $

In reliability and survival analysis, the survival function, also known as reliability function in the industrial context, is defined as $$S(t)= \Pr(T \gt t) = 1 - F(t),$$

S(t,a,b) := ''(1 - cdf_weibull(t, a, b));

$$S(t,a,b) := e^{-\left( \frac{t}{b} \right)^a}$$

Given a sample $(t_1, t_2, \ldots, t_n)$ with uncensored and right-censored data, the Kaplan-Meier estimator for the survival function is $$\hat{S}(t) = \prod_{t_j \leq t} \left( 1 - \frac{d_j}{n_j} \right),$$ where $d_j$ is the number of dying subjects at time $t_j$, and $n_j$ is the number of subjects still alive just before time $t_j$. An implementation of this algorithm follows,

/* Kaplan-Meier estimator */

load("descriptive")$
charfun2(z,l1,l2):= charfun(l1 <= z and z < l2)$

km(x) := block([s,ss,t,st,n,d,c,i:2,idx:1,hatS],
    s: sort(x),
    ss: length(x),
    t: append([0],
              listify(
                setify(
                  first(
                    transpose(
                      subsample(s,lambda([z], second(z)=1),1)))))),
    st: length(t),
    n: makelist(0,k,st),
    d: makelist(0,k,st),
    c: makelist(0,k,st),
    n[1]: ss,
    while i <= st do
        (   n[i]: n[i-1] - d[i-1] - c[i-1],
            while  idx <= ss and s[idx][1] <= t[i] do
                (   if s[idx][2] = 1
                        then d[i]: d[i] + 1
                    elseif s[idx][1] < t[i]
                        then n[i]: n[i] - 1
                        else c[i]: c[i] + 1,
                    idx: idx + 1),
            i: i + 1 ),
    hatS: 1- d / n,
    for k:2 thru st do
        hatS[k]: hatS[k-1] * hatS[k],
    sum(charfun2(z,t[i-1],t[i]) * hatS[i-1], i, 2, st)  + 
      charfun2(z,t[st],inf) * hatS[st] ) $

In order to see the Kaplan-Meier estimator in action, let's simulate a random sample of size 15 drawn from a Weibull distribution with parameters $a=10$ and $b=7$. Together with the time values we include a number to indicate whether the time is censored (0) or not (1). We consider observations longer than 8 to be right censored,

fpprintprec:4 $ /* for pretty printing */
(a: 10, b: 7, c: 8) $
x: random_weibull(a, b, 15) $
x: map(lambda([z], if z < c then [z,1] else [8,0]), x);

$$\left[ \left[ 7.578 , 1 \right] , \left[ 4.671 , 1 \right] , \left[ 7.163 , 1 \right] , \left[ 6.658 , 1 \right] , \left[ 8 , 0 \right] , \\ \left[ 6.905 , 1 \right] , \left[ 6.248 , 1 \right] , \left[ 5.441 , 1 \right] , \left[ 6.05 , 1 \right] , \left[ 6.045 , 1 \right] , \\ \left[ 7.158 , 1 \right] , \left[ 7.158 , 1 \right] , \left[ 7.566 , 1 \right] , \left[ 6.001 , 1 \right] , \left[ 5.537 , 1 \right] \right]$$

And here is the estimated survival step function,

hatS(z) := ''(km(x)) ;

$${{2\,{\it charfun}\left(7.566\leq z\land z<7.578\right)}\over{15}}+ {{{\it charfun}\left(7.163\leq z\land z<7.566\right)}\over{5}}+ \\{{4\, {\it charfun}\left(7.158\leq z\land z<7.163\right)}\over{15}}+ {{ {\it charfun}\left(7.158\leq z\land z<7.158\right)}\over{3}}+ \\ {{2\, {\it charfun}\left(6.905\leq z\land z<7.158\right)}\over{5}}+ {{7\, {\it charfun}\left(6.658\leq z\land z<6.905\right)}\over{15}}+ \\ {{8\, {\it charfun}\left(6.248\leq z\land z<6.658\right)}\over{15}}+ {{3\, {\it charfun}\left(6.05\leq z\land z<6.248\right)}\over{5}}+ \\ {{2\, {\it charfun}\left(6.045\leq z\land z<6.05\right)}\over{3}}+ {{11\, {\it charfun}\left(6.001\leq z\land z<6.045\right)}\over{15}}+ \\ {{4\, {\it charfun}\left(5.537\leq z\land z<6.001\right)}\over{5}}+ {{13\, {\it charfun}\left(5.441\leq z\land z<5.537\right)}\over{15}}+ \\ {{14\, {\it charfun}\left(4.671\leq z\land z<5.441\right)}\over{15}}+ {\it charfun}\left(0\leq z\land z<4.671\right)+ \\ {{{\it charfun}\left( 7.578\leq z\right)}\over{15}}$$

As an example, let's calculate the probability of a subject to be still alive at time $t=7$,

hatS(7);

$${{2}\over{5}}$$

And a plot of it,

draw2d(
  terminal         = png,
  grid             = true,
  background_color = orange,
  dimensions       = [400, 300],
  yrange           = [0, 1],
  title            = "Kaplan-Meier estimator",
  line_width       = 3,
  explicit(hatS(z), z, 0, 8),
  color            = white,
  line_width       = 1,
  explicit(S(z, a, b), z, 0, 8)) $

We have plotted in white color the true survival function, which we know, since our sample have been generated by simulation.

Let's now try to estimate parameters a and b from the sample. One of the standard methods is via the maximum likelihood. But in this case we'll use a method based on linear regression. See the following result,

kill(a,b)$
'S(z,a,b) = S(z,a,b);

$$S\left(z , a , b\right)=e^ {- \left({{z}\over{b}}\right)^{a} }$$

log(lhs(%)) = log(rhs(%));

$$\log S\left(z , a , b\right)=-\left({{z}\over{b}}\right)^{a}$$

-lhs(%) = -rhs(%);

$$-\log S\left(z , a , b\right)=\left({{z}\over{b}}\right)^{a}$$

log(lhs(%)) = ev(log(rhs(%)), logexpand = all);

$$\log \left(-\log S\left(z , a , b\right)\right)=a\,\left(\log z- \log b\right)$$

lhs(%) = expand(rhs(%));

$$\log \left(-\log S\left(z , a , b\right)\right)=a\,\log z-a\,\log b$$

According to this last result, there is a linear relationship between $\log z$ and $\log \left(-\log S\left(z , a , b\right)\right)$. We can estimate parameters $m=a$ and $n=-a\,\log b$ by simple linear regression and then obtain estimators for a and b.

We extract from the sample true life times,

x1: subsample(x, lambda([v], second(v)=1),1);

$$\pmatrix{7.578\cr 4.671\cr 7.163\cr 6.658\cr 6.905\cr 6.248\cr 5.441\cr 6.05\cr 6.045\cr 7.158\cr 7.158\cr 7.566\cr 6.001\cr 5.537 \cr }$$

Estimate the survival time for each observation,

y1: float(matrixmap(hatS,x1));

$$\pmatrix{0.06667\cr 0.9333\cr 0.2\cr 0.4667\cr 0.4\cr 0.5333\cr 0.8667\cr 0.6\cr 0.6667\cr 0.2667\cr 0.3333\cr 0.1333\cr 0.7333\cr 0.8\cr }$$

Build a table with the two lists, removing those pairs with zero survival probability to avoid problems with logaritms. This task is performed with function subsample defined in package descriptive,

w: subsample(addcol(x1,y1), lambda([v], second(v)#0.0),1,2);

$$\pmatrix{7.578&0.06667\cr 4.671&0.9333\cr 7.163&0.2\cr 6.658&0.4667 \cr 6.905&0.4\cr 6.248&0.5333\cr 5.441&0.8667\cr 6.05&0.6\cr 6.045& 0.6667\cr 7.158&0.2667\cr 7.158&0.3333\cr 7.566&0.1333\cr 6.001& 0.7333\cr 5.537&0.8\cr }$$

Now we transformate this matrix to obtain pairs of the form $\left( \log z, \log \left(-\log \hat{S}\left(z\right) \right) \right)$. Function transform_sample, also from package descriptive, is our friend here,

pts: transform_sample(w, [a,b], [log(a),log(-log(b))]) $

$$\pmatrix{2.025&0.9962\cr 1.541&-2.674\cr 1.969&0.4759\cr 1.896&- 0.2716\cr 1.932&-0.08742\cr 1.832&-0.4642\cr 1.694&-1.944\cr 1.8&- 0.6717\cr 1.799&-0.9027\cr 1.968&0.279\cr 1.968&0.09405\cr 2.024& 0.7006\cr 1.792&-1.171\cr 1.711&-1.5\cr }$$

We are now ready for a least square adjustment,

load (lsquares)$

le: float(lsquares_estimates(pts,
                             [u,v],
                             v=m*u+n,
                             [m,n],
                             initial=[3,3],
                             iprint=[-1,0]));

$$\left[ \left[ m=7.315 , n=-14.07 \right] \right]$$

Since $m=a$ and $n=-a\,\log b$, the estimator for a is

aa: rhs(le[1][1]) ;

$$7.315$$

And for b,

float(solve(−aa*log(b)= rhs(le[1][2])));

$$\left[ b=6.845 \right]$$

Remember that our sample was obtained by simulation with $a=10$ and $b=7$.

© 2017, TecnoStats.