Some algorithms to fit some reliability mixture models under - - PowerPoint PPT Presentation

Some algorithms to fit some reliability mixture models under censoring

Laurent Bordes Didier Chauveau University of Pau University of Orl´ eans COMPSTAT August 22-27, 2010

Laurent Bordes () Fitting some reliability mixture models 27 August 2010 1 / 35

Table of contents

1

Reliability mixture models

2

Some real data sets

3

Parametric EM-algorithm

4

Parametric stochastic EM-algorithm

5

Semiparametric stochastic EM-algorithm

Laurent Bordes () Fitting some reliability mixture models 27 August 2010 2 / 35

Reliability mixture models

Plan

1

Reliability mixture models

2

Some real data sets

3

Parametric EM-algorithm

4

Parametric stochastic EM-algorithm

5

Semiparametric stochastic EM-algorithm

Laurent Bordes () Fitting some reliability mixture models 27 August 2010 3 / 35

Reliability mixture models

About lifetimes

The lifetime data are assumed to come from a finite mixture of m component densities fj, j = 1, . . . , m, where fj(·) = f (·|ξj) ∈ F a parametric family indexed by a Euclidean parameter ξ. The lifetime density of an observation X may be written X ∼ g(x|θ) =

m

• j=1

λjf (x|ξj), where θ = (λ, ξ) = (λ1, . . . , λm, ξ1, . . . , ξm). Latent variable representation: X = YZ where Z ∼ Mult(1, λ) and (YZ|Z = j) ∼ f (·|ξj). For references on the broad literature of mixture models McLachlan and Peel (2000).

Laurent Bordes () Fitting some reliability mixture models 27 August 2010 4 / 35

Reliability mixture models

Right censored data

The censoring process is described by a random variable C with density function q, distribution function Q and survival function ¯

• Q. In the right

censoring setup the only available information is T = min(X, C), D = I(X ≤ C). The n lifetime data are x = (x1, . . . , xn) iid∼ g, associated to n censoring times c = (c1, . . . , cn) iid∼ C. The observations are thus (t, d) = ((t1, d1), . . . , (tn, dn)) , where ti = min(xi, ci) and di = I(xi ≤ ci).

Laurent Bordes () Fitting some reliability mixture models 27 August 2010 5 / 35

Reliability mixture models

Complete data choice

The observed data (t, d) depends on x which comes from a finite mixture ⇒ missing data are naturally associated to it. To these incomplete data are associated complete data which correspond to the situation where the component of origin zi ∈ {1, . . . , m} of each individual lifetime xi is known. The complete model at the level of (X, Z) is given by Pθ(Z = z) = λz and (X|Z = z) ∼ fz. With the right censoring process the complete data are (t, d, z), where z = (z1, . . . , zn). Remark. As in Chauveau (1995) the complete data can be (x, z) instead of (t, d, z).

Laurent Bordes () Fitting some reliability mixture models 27 August 2010 6 / 35

Reliability mixture models

Complete data pdf

Because we have: f c

θ (T = t, D = 1, Z = z)

= Pθ(Z = z) fθ(D = 1, T = t|Z = z) = λz fθ(C ≥ X, X = t|z) = λz Pθ(C ≥ t) fθ(X = t|z) = λz fz(t) ¯ Q(t), and similarly f c

θ (t, 0, z) = λz ¯

Fz(t)q(t), the complete data pdf is summarized by f c(t, d, z|θ) =

• λzf (t|ξz) ¯

Q(t) d λz ¯ F(t|ξz)q(t) 1−d .

Laurent Bordes () Fitting some reliability mixture models 27 August 2010 7 / 35

Some real data sets

Plan

1

Reliability mixture models

2

Some real data sets

3

Parametric EM-algorithm

4

Parametric stochastic EM-algorithm

5

Semiparametric stochastic EM-algorithm

Laurent Bordes () Fitting some reliability mixture models 27 August 2010 8 / 35

Some real data sets

Acute Myelogenous Leukemia survival data (Miller, 1997)

Group effect with two groups Sample size: 23 Censored lifetimes: 5

Variables Description time survival or censoring time status censoring status x maintenance chemotherapy given group scale estimation Maintained 63.3 Nonmaintained 25.1

50 100 150 1 2 3 4 5

Laurent Bordes () Fitting some reliability mixture models 27 August 2010 9 / 35

Some real data sets

Lifetimes of diesel engines fans (Nelson, 1982)

Time scale (1000s of hours) Sample size: 70 Censored lifetimes: 12

Variables Description time survival or censoring time status censoring status

2000 4000 6000 8000 10000 2 4 6 8 5000 10000 15000 0.0 0.2 0.4 0.6 0.8 1.0

Reliability estimation

time (1000 hours) R Kaplan−Meier 1−Weibull

Laurent Bordes () Fitting some reliability mixture models 27 August 2010 10 / 35

Parametric EM-algorithm

Plan

1

Reliability mixture models

2

Some real data sets

3

Parametric EM-algorithm

4

Parametric stochastic EM-algorithm

5

Semiparametric stochastic EM-algorithm

Laurent Bordes () Fitting some reliability mixture models 27 August 2010 11 / 35

Parametric EM-algorithm

Parametric EM-algorithm: complete data = (t, d, z)

Usual missing data framework (Dempster, Laird and Rubin, 1977) ⇒ define an EM algorithm that generates a sequence (θk)k=1,2,... (with arbitrary initial value θ0) by iteratively maximize Q(θ|θk) = E

• log f c(t, d, Z|θ) | t, d, θk

=

n

• i=1

E

• log f c(ti, di, Zi|θ) | ti, di, θk

. Calculation of Q(θ|θk) requires calculation of the following posterior probabilities pk

ij

:= P(Zi = j|ti, di, θk) = λk

j

• f (ti|ξk

j )

p

ℓ=1 λk ℓ f (ti|ξk ℓ )

di ¯ F(ti|ξk

j )

p

ℓ=1 λk ℓ ¯

F(ti|ξk

ℓ )

1−di . (1)

Laurent Bordes () Fitting some reliability mixture models 27 August 2010 12 / 35

Parametric EM-algorithm

Exponential lifetimes: complete data = (t, d, z)

EM algorithm: θk → θk+1

1 E-step:

Calculate the posterior probabilities pk

ij as in Equation (1),

for all i = 1, . . . , n and j = 1, . . . , m.

2 M-step:

Set λk+1

j

= 1 n

n

• i=1

pk

ij

for j = 1, . . . , m ξk+1

j

= n

i=1 pk ij di

n

i=1 pk ij ti

for j = 1, . . . , m.

Laurent Bordes () Fitting some reliability mixture models 27 August 2010 13 / 35

Parametric EM-algorithm

Simulation example

g(x) = λ1ξ1 exp(−ξ1x) + λ2ξ2 exp(−ξ2x) x > 0, with ξ1 = 1 − − −, ξ2 = 0.2 − − − and λ1 = 1/3 − − −.

200 400 600 800 1000 0.0 0.5 1.0 1.5

EM for RMM, n=200, 30% censored

iterations estimates rate 1 rate 2 lambda 1 50 100 150 200 0.0 0.5 1.0 1.5

EM for RMM, n=1000, 34.7% censored

iterations estimates rate 1 rate 2 lambda 1

Laurent Bordes () Fitting some reliability mixture models 27 August 2010 14 / 35

Parametric EM-algorithm

Application to AML data: be careful!

20 40 60 80 100 20 40 60 80

Scale (100 iterations)

Iterations 100 200 300 400 500 20 40 60 80

Scale (500 iterations)

Iterations 20 40 60 80 100 0.0 0.4 0.8

Lambda (100 iterations)

Iterations 100 200 300 400 500 0.0 0.4 0.8

Lambda (500 iterations)

Iterations Laurent Bordes () Fitting some reliability mixture models 27 August 2010 15 / 35

Parametric EM-algorithm

Parametric EM-algorithm: complete data = (x, z)

Complete data pdf f c(x, z) = λzfz(x). Then Q(θ|θk) = E

• log f c(X, Z|θ) | t, d, θk

=

n

• i=1

E

• log f c(Xi, Zi|θ) | ti, di, θk

. Calculation of Q(θ|θk) requires calculation of the following posterior pdf f k

i (x, j)

:= f (Xi = x, Zi = j|ti, di, θk) = λk

j

• I(x = ti)f (ti|ξk

j )

p

ℓ=1 λk ℓ f (ti|ξk ℓ )

di I(x > ti)f (x|ξk

j )

p

ℓ=1 λk ℓ ¯

F(ti|ξk

ℓ )

1−di .

Laurent Bordes () Fitting some reliability mixture models 27 August 2010 16 / 35

Parametric EM-algorithm

Exponential lifetimes: complete data = (x, z)

EM algorithm: θk → θk+1

1 E-step:

Calculate the posterior probabilities pk

ij as in Equation (1),

for all i = 1, . . . , n and j = 1, . . . , m.

2 M-step:

Set for j = 1, . . . , m λk+1

j

= 1 n

n

• i=1

pk

ij,

ξk+1

j

= n

i=1 pk ij

n

i=1

• ditipk

ij + (1 − di) λk

j (1+ξk j ti) exp(−ξk j i)

ξk

j

p

ℓ=1 λk ℓ exp(−ξk ℓ ti)

.

Laurent Bordes () Fitting some reliability mixture models 27 August 2010 17 / 35

Parametric EM-algorithm

Remarks about the parametric EM algorithms

+ Whatever the choice of complete data the M-step for the λjs always leads to explicit formula. − Q(θ|θk) depends strongly on the choice of the underlying parametric family F. − Except for exponential lifetimes, explicit maximizers of Q(θ|θk) are not reachable. − Maximizing Q(θ|θk) may be as complicated as maximizing the full likelihood function.

Laurent Bordes () Fitting some reliability mixture models 27 August 2010 18 / 35

Parametric stochastic EM-algorithm

Plan

1

Reliability mixture models

2

Some real data sets

3

Parametric EM-algorithm

4

Parametric stochastic EM-algorithm

5

Semiparametric stochastic EM-algorithm

Laurent Bordes () Fitting some reliability mixture models 27 August 2010 19 / 35

Parametric stochastic EM-algorithm

Parametric stochastic EM approach [1/2]

Idea by Celeux and Diebolt (1985, 1986): at each iteration add a stochastic step where the missing data are simulated according to their posterior probability distribution given the current value θk of the unknown parameter θ. What should be the complete data? It is enough to chose (t, d, z). pk

i = (pk i1, . . . , pk im) is the posterior probability vector associated to

• bservation i. Consider Z ∼ Mult(1, pk

i ) a multinomial distributed

random variable with parameters 1 and pk

i (i.e., Z ∈ {1, . . . m} with

probabilities P(Z = j) = pk

ij)).

Laurent Bordes () Fitting some reliability mixture models 27 August 2010 20 / 35

Parametric stochastic EM-algorithm

Parametric stochastic EM approach [2/2]

St-EM algorithm: θk → θk+1

1 E-step:

Calculate the posterior probabilities pk

ij as in Equation (1),

for all i = 1, . . . , n and j = 1, . . . , m.

2 Stochastic step:

Simulate Z k

i ∼ Mult(1, pk i ), i = 1, . . . , n, and

define the subsets χk

j = {i ∈ {1, . . . , n} : Z k i = j},

j = 1, . . . , m. (2)

3 M-step:

For each component j ∈ {1, . . . , m} λk+1

j

= Card(χk

j )/n,

and ξk+1

j

= arg max

ξ∈Ξ Lj(ξ),

where Lj(ξ) =

• i∈χk

j

(f (ti|ξ))di(¯ F(ti|ξ))1−di.

Laurent Bordes () Fitting some reliability mixture models 27 August 2010 21 / 35

Parametric stochastic EM-algorithm

Exponential mixture example (n = 200)

g(x) = λξ1 exp(−ξ1x) + (1 − λ)ξ2 exp(−ξ2x) x > 0, with λ = 1/3, ξ1 = 1 and ξ2 = 1/5.

20 40 60 80 100 0.1 0.2 0.3 0.4 0.5

lambda_1 estimation

iterations 20 40 60 80 100 0.0 0.5 1.0 1.5 2.0

xi's estimation

iterations

Laurent Bordes () Fitting some reliability mixture models 27 August 2010 22 / 35

Parametric stochastic EM-algorithm

Application to engine fans (two exponentials)

100 300 500 0.0 0.2 0.4

lambda_1 estimation

iterations 100 300 500 0e+00 2e−04 4e−04

xi's estimation

iterations

Laurent Bordes () Fitting some reliability mixture models 27 August 2010 23 / 35

Parametric stochastic EM-algorithm

Application to engine fans (two Weibulls)

100 200 300 400 500 0.0 0.4 0.8

lambda estimation

iterations 100 200 300 400 500 4000 8000

Scales estimation

iterations 100 200 300 400 500 1 2 3 4 5 6

Shapes estimation

iterations 5000 10000 15000 0.0 0.4 0.8

Empirical vs RMM reliab.

time (1000 hours) R Laurent Bordes () Fitting some reliability mixture models 27 August 2010 24 / 35

Semiparametric stochastic EM-algorithm

Plan

1

Reliability mixture models

2

Some real data sets

3

Parametric EM-algorithm

4

Parametric stochastic EM-algorithm

5

Semiparametric stochastic EM-algorithm

Laurent Bordes () Fitting some reliability mixture models 27 August 2010 25 / 35

Semiparametric stochastic EM-algorithm

A semiparametric reliability mixture model (SRMM)

g(x|θ) = λ1f (x) + λ2ξf (ξx) x > 0, where θ = (λ, ξ, f ) ∈ (0, 1) × R+

∗ × F; F is a family of pdf.

Interpretation: accelerated lifetime model for grouped data with two groups and unobserved group label. Example: λ1 = 0.3, ξ = 0.1 and f ∼ LN(1, 0.5).

10 20 30 40 50 60 0.00 0.05 0.10 0.15

Laurent Bordes () Fitting some reliability mixture models 27 August 2010 26 / 35

Semiparametric stochastic EM-algorithm

Identifiabillity of θ

Question: how to chose F to obtain

• ∀x > 0

g(x|θ) = g(x|θ′)

• ⇒ θ = θ′?

Hard question. . . partial answer in Bordes, Mottelet and Vandekerkhove (2006) and in Hunter, Wang and Hettmansperger (2007): if F is a subset

• f pdf f such that x → exf (ex) is symmetric then identifiability holds!

Laurent Bordes () Fitting some reliability mixture models 27 August 2010 27 / 35

Semiparametric stochastic EM-algorithm

Stochastic EM algorithm for the SRMM [1/4]

Notations: f is the unknown pdf, write ¯ F the reliability function and α = f /¯ F the failure rate. From (t, d): ¯ F is nonparametrically estimated by the Kaplan-Meier estimator, α is nonparametrically estimated by smoothing the Nelson-Aalen estimator. Because λk, ξk, ¯ F k and αk are estimates of λ, ξ, ¯ F and α at step k we have: pk

ij

:= P(Zi = j|ti, di, θk) =

• αk(ti)¯

F k(ti) p

ℓ=1 λk ℓ αk(ti)¯

F k(ti) di λk

j ¯

F k(ti) p

ℓ=1 λk ℓ ¯

F k(ti) 1−di , where the pdf f is estimated by f k = αk ¯ F k.

Laurent Bordes () Fitting some reliability mixture models 27 August 2010 28 / 35

Semiparametric stochastic EM-algorithm

Stochastic EM algorithm for the SRMM [2/4]

1 Posterior probabilities calculation: for each item i ∈ {1, . . . , n}:

if di = 0 then pk

i1 =

λk ¯ F k(ti) λk ¯ F k(ti) + (1 − λk)¯ F k(ξkti), else pk

i1 =

λkαk(ti)¯ F k(ti) λkαk(ti)¯ F k(ti) + (1 − λk)ξkαk(ξkti)¯ F k(ξkti). Set pk

i = (pk i1, 1 − pk i1).

2 Stochastic step: for each item i ∈ {1, . . . , n} simulate

Z k

i ∼ Mult(1, pk i ). Then define subsets

χk

j = {i ∈ {1, . . . , n}; Z k i = j}

for j = 1, 2.

Laurent Bordes () Fitting some reliability mixture models 27 August 2010 29 / 35

Semiparametric stochastic EM-algorithm

Stochastic EM algorithm for the SRMM [3/4]

Facts: ξ = E(X|Z=1)

E(X|Z=2) and if Sj(s) = P(X > s|Z = j) then

E(X|Z = j) = +∞ Sj(s)ds.

3 Update the euclidean parameters λ and ξ:

λk+1 = Card(χk

1)

n , ξk+1 = τ k

1

ˆ Sk

1 (s)ds

τ k

2

ˆ Sk

2 (s)ds

, where ˆ Sk

j is the Kaplan-Meier estimator for the subpopulation

{(tℓ, dℓ); ℓ ∈ χk

j } and τ k j = maxℓ∈χk

j tℓ. Laurent Bordes () Fitting some reliability mixture models 27 August 2010 30 / 35

Semiparametric stochastic EM-algorithm

Stochastic EM algorithm for the SRMM [4/4]

Fact: if X comes from component two (i.e. if Z = 2), then ξX ∼ f .

4 Update the functional parameters α and ¯

F: set tk = (tk

1 , . . . , tk n ) be

the order statistic from {ti; i ∈ χk

1} ∪ {ξkti; i ∈ χk 2}; write

dk = (dk

1 , . . . , dk n ) the corresponding censoring indicators.

αk+1(s) =

n

• i=1

1 hK s − tk

i

h

• dk

i

n − i + 1, ¯ F k+1(s) =

• i:tk

i ≤s

• 1 −

dk

i

n − i + 1

• ,

where K is a kernel function and h a bandwidth. Remark: in practice the choice of both K and h is important!

Laurent Bordes () Fitting some reliability mixture models 27 August 2010 31 / 35

Semiparametric stochastic EM-algorithm

Example: g(x) = 0.3f (x) + 0.7ξf (ξx), f ∼ LN(1, 0.5). Simulated sample: n = 100 with 0% of censoring.

20 40 60 80 100 0.1 0.2 0.3 0.4 0.5 Estimation of the mixture proportion Iterations lambda 20 40 60 80 100 0.00 0.05 0.10 0.15 0.20 Estimation of the scale parameter Iterations xi 2 4 6 8 10 0.0 0.2 0.4 0.6 0.8 1.0 Density estimation time density

Laurent Bordes () Fitting some reliability mixture models 27 August 2010 32 / 35

Semiparametric stochastic EM-algorithm

Example (continued): n = 200 and 10% of censoring.

2 4 6 8 10 12 0.0 0.1 0.2 0.3 0.4 0.5

scaling

iterations 2 4 6 8 10 12 0.0 0.1 0.2 0.3 0.4 0.5

weight of component 1

iterations 10 20 30 40 50 60 0.0 0.2 0.4 0.6 0.8 1.0

Survival ft

time f(x) estimate true 10 20 30 40 50 60 0.00 0.01 0.02 0.03 0.04

Density

time density estimate true

Laurent Bordes () Fitting some reliability mixture models 27 August 2010 33 / 35

Semiparametric stochastic EM-algorithm

Conclusions

All the algorithms introduced here have been/will be implemented in the publicly available package mixtools by Benaglia, Chauveau, Hunter and Young (2009) for the R statistical software (R Development Core Team, 2009). Asymptotic variances of the parametric St-EM estimators can be derived following Nielsen (2000). Many tuning parameters to improve. As an example, a local bandwidth choice should improve the semiparametric St-EM algorithm.

Laurent Bordes () Fitting some reliability mixture models 27 August 2010 34 / 35

Semiparametric stochastic EM-algorithm

Thanks!

Laurent Bordes () Fitting some reliability mixture models 27 August 2010 35 / 35