Package {PScr}

Maximum likelihood estimation based on EM algorithm for the Power Series cure rate model

Description

This function provides the maximum likelihood estimation based on the EM algorithm for the Power Series cure rate model

Usage

EM.PScr(t, delta, z, model = 1, dist = 1, max.iter = 1000, 
           prec = 1e-04)

Arguments

Details

The EM algorithm for the model is implemented as in Gallardo et al. (2017).

Value

Author(s)

Diego I. Gallardo and Reza Azimi

References

Azimi, R, Esmailian, M, Gallardo DI and Gomez HJ. (2022). A New Cure Rate Model Based on Flory-Schulz Distribution: Application to the Cancer Data. Mathematics 10, 4643

Gallardo DI, Gomez YM and De Castro M. (2018). A flexible cure rate model based on the polylogarithm distribution. Journal of Statistical Computation and Simulation 88 (11), 2137-2149

Gallardo DI, Gomez YM, Gomez HJ, Gallardo-Nelson MJ, Bourguignon M. (2022) The slash half-normal distribution applied to a cure rate model with application to bone marrow transplantation. Mathematics, Submitted.

Gallardo DI, Romeo JS and Meyer R. (2017). A simplified estimation procedure based on the EM algorithm for the power series cure rate model. Communications in Statistics-Simulation and Computation 46 (8), 6342-6359.

Examples

require(mstate)
data(ebmt4)
attach(ebmt4)
t = srv / 365.25 # Time in years
delta=srv.s
prophy=as.factor(proph)
year2=ifelse(year=="1985-1989",0,1)
z=t(model.matrix(~proph-1))
#Computes the estimation for Poisson-Slash half-normal cure rate model
EM.PScr(t, delta, z, model=1, dist=1)
#Computes the estimation for Flory-Schulz-Slash half-normal cure rate model
EM.PScr(t, delta, z, model=6, dist=1)

Simulation of Right-Censored Data from a Power Series Cure Rate Model

Description

This function generates right-censored survival data from a power series cure rate model. The number of latent competing causes follows a power series distribution (Poisson, Bernoulli, negative binomial, geometric or polylogarithm) and the time-to-event of each cause follows a Weibull, gamma or Birnbaum-Saunders distribution. The generated data can be analysed with EM.PScr.

Usage

rcura(n, beta, z = NULL, q = NULL, alpha, sigma, C, model, dist)

Arguments

Details

Let M_i be the number of latent causes of individual i, i = 1, \ldots, n, with power parameter \theta_i = g(z_i^{\top}\beta). Given M_i > 0, the time-to-event is T_i = \min(W_{i1}, \ldots, W_{iM_i}), where the W_{ij} are independent and identically distributed with the distribution selected in dist. Individuals with M_i = 0 are cured and have T_i = \infty. The observed time is t_i = \min(T_i, C) and the failure indicator is \delta_i = I(T_i \leq C).

The distributions available for the latent causes, the link functions and the corresponding cure fractions p_0 = P(M = 0) are:

The time-to-event distributions are parameterized by a shape alpha and a scale sigma: "weibull" uses rweibull, "gamma" uses rgamma and "bs" uses rbisa. This is the same parameterization used by EM.PScr, so the simulated data can be fitted with the matching dist code.

The cure fraction and the percentage of censoring of a simulated sample depend on beta, alpha, sigma and C. If the time-to-event distribution has a heavy tail relative to C, a fraction of the susceptible individuals is censored before the event occurs, so the Kaplan-Meier curve does not level off at p_0 and the cure fraction is difficult to identify. The empirical cure fraction and the proportion of censored observations of each sample are returned as attributes.

Value

A data.frame with n rows and two columns:

The data frame has the following attributes, which can be extracted with attr:

Note

The function uses the random number generator, so results are reproducible only after calling set.seed. The computing time grows linearly with n.

Correspondence with the codes of EM.PScr: model = "poisson" corresponds to model = 1, "binneg" and "geometrica" to model = 3 (with q = 1 in the geometric case), "bernoulli" to model = 4 and "polilog" to model = 5; dist = "weibull", "gamma" and "bs" correspond to dist = 2, 3 and 4, respectively. The logarithmic and Flory-Schulz models and the slash half-normal and log-normal time-to-event distributions available in EM.PScr are not implemented in rcura.

Author(s)

Daniel Jana, Yolanda Gomez and Diego Gallardo

References

Birnbaum ZW and Saunders SC. (1969). A new family of life distributions. Journal of Applied Probability 6 (2), 319-327.

Cancho VG, Louzada F and Ortega EMM. (2013). The power series cure rate model: an application to a cutaneous melanoma data. Communications in Statistics - Simulation and Computation 42 (3), 586-602. doi:10.1080/03610918.2011.639971

Gallardo DI, Gomez YM and De Castro M. (2018). A flexible cure rate model based on the polylogarithm distribution. Journal of Statistical Computation and Simulation 88 (11), 2137-2149.

Gallardo DI, Romeo JS and Meyer R. (2017). A simplified estimation procedure based on the EM algorithm for the power series cure rate model. Communications in Statistics - Simulation and Computation 46 (8), 6342-6359. doi:10.1080/03610918.2016.1202276

Yakovlev AY and Tsodikov AD. (1996). Stochastic Models of Tumor Latency and Their Biostatistical Applications. World Scientific, Singapore.

See Also

EM.PScr

Examples

set.seed(2026)

## Poisson cure rate model with Weibull times and no covariates.
## With theta = -log(0.7) the cure fraction is p0 = exp(-theta) = 0.7
theta <- -log(0.7)
dat <- rcura(n = 500, beta = log(theta), alpha = 2, sigma = 3,
             C = 8, model = "poisson", dist = "weibull")
head(dat[dat$delta == 1, ])  # first observed events
attr(dat, "p0_empirico")     # empirical cure fraction (approximately 0.7)
attr(dat, "cens_pct")        # proportion of censored observations

## The Kaplan-Meier curve levels off near the cure fraction
fit <- survival::survfit(survival::Surv(t, delta) ~ 1, data = dat)
plot(fit, conf.int = FALSE, xlab = "Time", ylab = "Survival probability")
abline(h = 0.7, lty = 2)

## Negative binomial model (q = 2) with gamma times and cure fraction 0.4.
## Since p0 = (1 - theta)^q, we have theta = 1 - p0^(1 / q)
q <- 2
theta <- 1 - 0.4^(1 / q)
dat2 <- rcura(n = 200, beta = qlogis(theta), q = q, alpha = 2, sigma = 1.4,
              C = 8, model = "binneg", dist = "gamma")
mean(dat2$delta == 0)

## Intercept and a binary covariate (covariates in rows, individuals in columns)
n <- 200
z <- rbind(1, rbinom(n, 1, 0.5))
dat3 <- rcura(n = n, beta = c(-0.5, 0.7), z = z, alpha = 0.5, sigma = 2.5,
              C = 8, model = "poisson", dist = "bs")
table(dat3$delta)


## Fit the first simulated data set with EM.PScr (model = 1: Poisson, dist = 2: Weibull)
fit.em <- EM.PScr(t = dat$t, delta = dat$delta,
                  z = matrix(1, nrow = 1, ncol = nrow(dat)), model = 1, dist = 2)
fit.em$estimate

## Polylogarithm model
dat4 <- rcura(n = 200, beta = 0.5, q = 2, alpha = 2, sigma = 3,
              C = 8, model = "polilog", dist = "weibull")
attr(dat4, "p0_empirico")