St-EM algorithm for Reliability Mixture Models (RMM) of Weibull with right Censoring
weibullRMM_SEM.RdParametric Stochastic EM (St-EM) algorithm for univariate finite mixture of Weibull distributions with randomly right censored data.
Usage
weibullRMM_SEM(x, d = NULL, lambda = NULL, shape = NULL, scale = NULL,
k = 2, maxit = 200, maxit.survreg = 200, epsilon = 1e-03,
averaged = TRUE, verb = FALSE)Arguments
- x
A vector of \(n\) real positive lifetime (possibly censored) durations. If
dis notNULLthen a vector of random censoring timescoccured, so that \(x= min(x,c)\) and \(d = I(x <= c)\).- d
The vector of censoring indication, where 1 means observed lifetime data, and 0 means censored lifetime data.
- lambda
Initial value of mixing proportions. If
NULL, thenlambdais set torep(1/k,k).- shape
Initial value of Weibull component shapes, all set to 1 if
NULL.- scale
Initial value of Weibull component scales, all set to 1 if
NULL.- k
Number of components of the mixture.
- maxit
The number of iterations allowed, since for St-EM algorithms convergence is not based on stabilization, exactly
maxititerations are performed (see Bordes L. and Chauveau D. (2016) reference below).- maxit.survreg
The number of iterations allowed in the computations of the MLE for censored weibull data from the
survivalpackage (see Bordes L. and Chauveau D. (2016) reference below).- epsilon
Tolerance parameter used in the numerical computations of the MLE for censored weibull data by
survregfrom thesurvivalpackage (see Bordes L. and Chauveau D. (2016) reference below).- averaged
The way of updating parameters at each iteration: if
TRUE, current values of the parameters are obtained by averaging the sequence (see Bordes L. and Chauveau D. (2016) reference below).- verb
If TRUE, print updates for every iteration of the algorithm as it runs
Details
This St-EM algorithm calls functions from the survival package to compute
parametric MLE for censored weibull data.
Value
weibullRMM_SEM returns a list of class "mixEM" with the following items:
- x
The input data.
- d
The input censoring indicator.
- lambda
The estimates for the mixing proportions.
- scale
The estimates for the Weibull component scales.
- shape
The estimates for the Weibull component shapes.
- loglik
The log-likelihood value at convergence of the algorithm.
- posterior
An \(n\times k\) matrix of posterior probabilities for observation, after convergence of the algorithm.
- all.loglik
The sequence of log-likelihoods over iterations.
- all.lambda
The sequence of mixing proportions over iterations.
- all.scale
The sequence of component scales over iterations.
- all.shape
The sequence of component shapes over iterations.
- ft
A character vector giving the name of the function called.
See also
Related functions:
plotweibullRMM, summary.mixEM.
Other models and algorithms for censored lifetime data
(name convention is model_algorithm):
expRMM_EM,
spRMM_SEM.
References
Bordes, L., and Chauveau, D. (2016), Stochastic EM algorithms for parametric and semiparametric mixture models for right-censored lifetime data, Computational Statistics, Volume 31, Issue 4, pages 1513-1538. https://link.springer.com/article/10.1007/s00180-016-0661-7
Examples
n = 500 # sample size
m = 2 # nb components
lambda=c(0.4, 0.6)
shape <- c(0.5,5); scale <- c(1,20) # model parameters
set.seed(321)
x <- rweibullmix(n, lambda, shape, scale) # iid ~ weibull mixture
cs=runif(n,0,max(x)+10) # iid censoring times
t <- apply(cbind(x,cs),1,min) # censored observations
d <- 1*(x <= cs) # censoring indicator
## set arbitrary or "reasonable" (e.g., data-driven) initial values
l0 <- rep(1/m,m); sh0 <- c(1, 2); sc0 <- c(2,10)
# Stochastic EM algorithm
a <- weibullRMM_SEM(t, d, lambda = l0, shape = sh0, scale = sc0, maxit = 200)
#> number of iterations = 200
summary(a) # Parameters estimates etc
#> summary of weibullRMM_SEM object:
#> comp 1 comp 2
#> lambda 0.396100 0.60390
#> shape 0.497675 5.18371
#> scale 0.620329 20.20621
#> loglik at estimate: -960.1892
#> 25 % of the data right censored
plotweibullRMM(a) # plot of St-EM sequences
plot(a, which=2) # or equivalently, S3 method for "mixEM" object