# St-EM algorithm for Reliability Mixture Models (RMM) of Weibull with right Censoring Parametric Stochastic EM (St-EM) algorithm for univariate finite mixture of Weibull distributions with randomly right censored data. ```r 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 `d` is not `NULL` then a vector of random censoring times `c` occured, 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`, then `lambda` is set to `rep(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 `maxit` iterations 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 `survival` package (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 `survreg` from the `survival` package (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. ## Returns `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 x 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](https://link.springer.com/article/10.1007/s00180-016-0661-7) ## Author(s) Didier Chauveau ## Examples ```r 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) summary(a) # Parameters estimates etc plotweibullRMM(a) # plot of St-EM sequences plot(a, which=2) # or equivalently, S3 method for "mixEM" object ```

weibullRMM_SEM function