| Title: | Split-Population (Cure / Mover-Stayer) Survival Models |
|---|---|
| Description: | Maximum-likelihood estimation of split-population (cure / mover-stayer) survival models: an accelerated failure-time regression for event timing among "movers", combined with a logistic regression on the probability of belonging to the immune "stayer" population. Five baseline timing distributions are provided -- log-logistic, Weibull, log-normal, gamma, and the generalized gamma that nests the other four -- following Schmidt & Witte (1989, Journal of Econometrics) and Yamaguchi (1992, 1998, Sociological Methodology). This is an R translation of a set of Stata ml programs, with the log-likelihood corrected to match the published model and verified by simulation against known parameters. |
| Authors: | Nobutaka Fukuda [aut, cre] |
| Maintainer: | Nobutaka Fukuda <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 0.1.0 |
| Built: | 2026-10-06 11:39:22 UTC |
| Source: | https://github.com/nobifukuda/splitpopsurv |
Maximum-likelihood estimation of a split-population (cure) survival model with a gamma baseline for the "mover" population. R translation of the Stata program 'SphGam'; see 'docs/manual.html' for the model and a note on a likelihood correction relative to the original Stata code.
fit_splitpop_gamma( hform, pform, data, time, event, group, start = NULL, method = "NR" )fit_splitpop_gamma( hform, pform, data, time, event, group, start = NULL, method = "NR" )
hform |
One-sided formula (e.g. '~ x1 + x2') for the H_regression equation: covariates for the log-logistic timing distribution. |
pform |
One-sided formula for the P_regression equation: covariates for the logit on the stayer ("cure") probability. |
data |
A 'data.frame' containing 'time', 'event', 'group', and every variable referenced in 'hform'/'pform'. |
time |
Character: name of the duration column ($ML_y1 in Stata). |
event |
Character: name of the 0/1 event column, 1 = failure observed, 0 = censored ($ML_y2 in Stata). |
group |
Character: name of a 0/1 column that flips which side of the cure-probability logit is used ($ML_y3 in Stata) – see the manual. If there is no such distinction in your data, pass a column of all 1s. |
start |
Optional numeric starting vector; if 'NULL', a default is constructed automatically. |
method |
Optimizer passed to [maxLik::maxLik()]; default '"NR"' (Newton-Raphson). Try '"BFGS"' if that fails to converge. |
Note: unlike the other four models, the H_regression linear index here is used directly as a rate (no 'exp()' transform), so it must stay positive; the default 'start' sets a positive intercept and zero slopes to keep it safe regardless of covariate values.
A 'maxLik' object; use 'summary()', 'coef()', 'logLik()', 'vcov()'.
set.seed(1) n <- 150 x1 <- rnorm(n); x2 <- rbinom(n, 1, 0.5) cured <- rbinom(n, 1, plogis(0.2 + 0.5 * x2)) t_latent <- (-log(runif(n)) / exp(0.5 - 0.4 * x1))^(1 / 1.3) censor_time <- rexp(n, rate = 0.2) time <- pmax(ifelse(cured == 1, censor_time, pmin(t_latent, censor_time)), 1e-3) event <- ifelse(cured == 1, 0, as.numeric(t_latent <= censor_time)) mydata <- data.frame(time = time, event = event, x1 = x1, x2 = x2, group = 1) fit <- fit_splitpop_gamma(~x1, ~x2, mydata, "time", "event", "group", method = "BFGS") summary(fit)set.seed(1) n <- 150 x1 <- rnorm(n); x2 <- rbinom(n, 1, 0.5) cured <- rbinom(n, 1, plogis(0.2 + 0.5 * x2)) t_latent <- (-log(runif(n)) / exp(0.5 - 0.4 * x1))^(1 / 1.3) censor_time <- rexp(n, rate = 0.2) time <- pmax(ifelse(cured == 1, censor_time, pmin(t_latent, censor_time)), 1e-3) event <- ifelse(cured == 1, 0, as.numeric(t_latent <= censor_time)) mydata <- data.frame(time = time, event = event, x1 = x1, x2 = x2, group = 1) fit <- fit_splitpop_gamma(~x1, ~x2, mydata, "time", "event", "group", method = "BFGS") summary(fit)
Maximum-likelihood estimation of a split-population (cure) survival model with a generalized gamma baseline (Prentice 1974 / Yamaguchi & Ferguson 1995, note 10) for the "mover" population – the family that nests the Weibull (kappa=1), log-normal (kappa=0), and gamma (sigma=1) models above. R translation of the Stata program 'SphGGam' ('d0' method); see 'docs/manual.html' for the model and a note on a likelihood correction relative to the original Stata code.
fit_splitpop_ggamma( hform, pform, data, time, event, group, start = NULL, method = "NR" )fit_splitpop_ggamma( hform, pform, data, time, event, group, start = NULL, method = "NR" )
hform |
One-sided formula (e.g. '~ x1 + x2') for the H_regression equation: covariates for the log-logistic timing distribution. |
pform |
One-sided formula for the P_regression equation: covariates for the logit on the stayer ("cure") probability. |
data |
A 'data.frame' containing 'time', 'event', 'group', and every variable referenced in 'hform'/'pform'. |
time |
Character: name of the duration column ($ML_y1 in Stata). |
event |
Character: name of the 0/1 event column, 1 = failure observed, 0 = censored ($ML_y2 in Stata). |
group |
Character: name of a 0/1 column that flips which side of the cure-probability logit is used ($ML_y3 in Stata) – see the manual. If there is no such distinction in your data, pass a column of all 1s. |
start |
Optional numeric starting vector; if 'NULL', a default is constructed automatically. |
method |
Optimizer passed to [maxLik::maxLik()]; default '"NR"' (Newton-Raphson). Try '"BFGS"' if that fails to converge. |
Note: '|kappa| < 0.01' switches to a log-normal-equivalent branch, a genuine kink in the likelihood surface (preserved from the Stata guard). If 'summary()' reports a non-finite standard error, the optimizer likely landed near that threshold – try a different starting 'kappa'.
A 'maxLik' object; use 'summary()', 'coef()', 'logLik()', 'vcov()'.
set.seed(1) n <- 150 x1 <- rnorm(n); x2 <- rbinom(n, 1, 0.5) cured <- rbinom(n, 1, plogis(0.2 + 0.5 * x2)) t_latent <- (-log(runif(n)) / exp(0.5 - 0.4 * x1))^(1 / 1.3) censor_time <- rexp(n, rate = 0.2) time <- pmax(ifelse(cured == 1, censor_time, pmin(t_latent, censor_time)), 1e-3) event <- ifelse(cured == 1, 0, as.numeric(t_latent <= censor_time)) mydata <- data.frame(time = time, event = event, x1 = x1, x2 = x2, group = 1) fit <- fit_splitpop_ggamma(~x1, ~x2, mydata, "time", "event", "group", method = "BFGS") summary(fit)set.seed(1) n <- 150 x1 <- rnorm(n); x2 <- rbinom(n, 1, 0.5) cured <- rbinom(n, 1, plogis(0.2 + 0.5 * x2)) t_latent <- (-log(runif(n)) / exp(0.5 - 0.4 * x1))^(1 / 1.3) censor_time <- rexp(n, rate = 0.2) time <- pmax(ifelse(cured == 1, censor_time, pmin(t_latent, censor_time)), 1e-3) event <- ifelse(cured == 1, 0, as.numeric(t_latent <= censor_time)) mydata <- data.frame(time = time, event = event, x1 = x1, x2 = x2, group = 1) fit <- fit_splitpop_ggamma(~x1, ~x2, mydata, "time", "event", "group", method = "BFGS") summary(fit)
Maximum-likelihood estimation of a split-population (cure) survival model with a log-logistic baseline for the "mover" (susceptible) population and a logistic regression on the probability of being a "stayer" (immune to the event). R translation of the Stata program 'SphLog'; see 'docs/manual.html' for the model, formulas, and a note on a likelihood correction relative to the original Stata code.
fit_splitpop_loglogistic( hform, pform, data, time, event, group, start = NULL, method = "NR" )fit_splitpop_loglogistic( hform, pform, data, time, event, group, start = NULL, method = "NR" )
hform |
One-sided formula (e.g. '~ x1 + x2') for the H_regression equation: covariates for the log-logistic timing distribution. |
pform |
One-sided formula for the P_regression equation: covariates for the logit on the stayer ("cure") probability. |
data |
A 'data.frame' containing 'time', 'event', 'group', and every variable referenced in 'hform'/'pform'. |
time |
Character: name of the duration column ($ML_y1 in Stata). |
event |
Character: name of the 0/1 event column, 1 = failure observed, 0 = censored ($ML_y2 in Stata). |
group |
Character: name of a 0/1 column that flips which side of the cure-probability logit is used ($ML_y3 in Stata) – see the manual. If there is no such distinction in your data, pass a column of all 1s. |
start |
Optional numeric starting vector; if 'NULL', a default is constructed automatically. |
method |
Optimizer passed to [maxLik::maxLik()]; default '"NR"' (Newton-Raphson). Try '"BFGS"' if that fails to converge. |
A 'maxLik' object; use 'summary()', 'coef()', 'logLik()', 'vcov()'.
set.seed(1) n <- 150 x1 <- rnorm(n); x2 <- rbinom(n, 1, 0.5) cured <- rbinom(n, 1, plogis(0.2 + 0.5 * x2)) t_latent <- (-log(runif(n)) / exp(0.5 - 0.4 * x1))^(1 / 1.3) censor_time <- rexp(n, rate = 0.2) time <- pmax(ifelse(cured == 1, censor_time, pmin(t_latent, censor_time)), 1e-3) event <- ifelse(cured == 1, 0, as.numeric(t_latent <= censor_time)) mydata <- data.frame(time = time, event = event, x1 = x1, x2 = x2, group = 1) fit <- fit_splitpop_loglogistic(~x1, ~x2, mydata, "time", "event", "group", method = "BFGS") summary(fit)set.seed(1) n <- 150 x1 <- rnorm(n); x2 <- rbinom(n, 1, 0.5) cured <- rbinom(n, 1, plogis(0.2 + 0.5 * x2)) t_latent <- (-log(runif(n)) / exp(0.5 - 0.4 * x1))^(1 / 1.3) censor_time <- rexp(n, rate = 0.2) time <- pmax(ifelse(cured == 1, censor_time, pmin(t_latent, censor_time)), 1e-3) event <- ifelse(cured == 1, 0, as.numeric(t_latent <= censor_time)) mydata <- data.frame(time = time, event = event, x1 = x1, x2 = x2, group = 1) fit <- fit_splitpop_loglogistic(~x1, ~x2, mydata, "time", "event", "group", method = "BFGS") summary(fit)
Maximum-likelihood estimation of a split-population (cure) survival model with a log-normal baseline for the "mover" population. R translation of the Stata program 'SphNom'; see 'docs/manual.html' for the model and a note on a likelihood correction relative to the original Stata code.
fit_splitpop_lognormal( hform, pform, data, time, event, group, start = NULL, method = "NR" )fit_splitpop_lognormal( hform, pform, data, time, event, group, start = NULL, method = "NR" )
hform |
One-sided formula (e.g. '~ x1 + x2') for the H_regression equation: covariates for the log-logistic timing distribution. |
pform |
One-sided formula for the P_regression equation: covariates for the logit on the stayer ("cure") probability. |
data |
A 'data.frame' containing 'time', 'event', 'group', and every variable referenced in 'hform'/'pform'. |
time |
Character: name of the duration column ($ML_y1 in Stata). |
event |
Character: name of the 0/1 event column, 1 = failure observed, 0 = censored ($ML_y2 in Stata). |
group |
Character: name of a 0/1 column that flips which side of the cure-probability logit is used ($ML_y3 in Stata) – see the manual. If there is no such distinction in your data, pass a column of all 1s. |
start |
Optional numeric starting vector; if 'NULL', a default is constructed automatically. |
method |
Optimizer passed to [maxLik::maxLik()]; default '"NR"' (Newton-Raphson). Try '"BFGS"' if that fails to converge. |
A 'maxLik' object; use 'summary()', 'coef()', 'logLik()', 'vcov()'.
set.seed(1) n <- 150 x1 <- rnorm(n); x2 <- rbinom(n, 1, 0.5) cured <- rbinom(n, 1, plogis(0.2 + 0.5 * x2)) t_latent <- (-log(runif(n)) / exp(0.5 - 0.4 * x1))^(1 / 1.3) censor_time <- rexp(n, rate = 0.2) time <- pmax(ifelse(cured == 1, censor_time, pmin(t_latent, censor_time)), 1e-3) event <- ifelse(cured == 1, 0, as.numeric(t_latent <= censor_time)) mydata <- data.frame(time = time, event = event, x1 = x1, x2 = x2, group = 1) fit <- fit_splitpop_lognormal(~x1, ~x2, mydata, "time", "event", "group", method = "BFGS") summary(fit)set.seed(1) n <- 150 x1 <- rnorm(n); x2 <- rbinom(n, 1, 0.5) cured <- rbinom(n, 1, plogis(0.2 + 0.5 * x2)) t_latent <- (-log(runif(n)) / exp(0.5 - 0.4 * x1))^(1 / 1.3) censor_time <- rexp(n, rate = 0.2) time <- pmax(ifelse(cured == 1, censor_time, pmin(t_latent, censor_time)), 1e-3) event <- ifelse(cured == 1, 0, as.numeric(t_latent <= censor_time)) mydata <- data.frame(time = time, event = event, x1 = x1, x2 = x2, group = 1) fit <- fit_splitpop_lognormal(~x1, ~x2, mydata, "time", "event", "group", method = "BFGS") summary(fit)
Maximum-likelihood estimation of a split-population (cure) survival model with a Weibull baseline for the "mover" population. R translation of the Stata program 'SphWieb'; see 'docs/manual.html' for the model and a note on a likelihood correction relative to the original Stata code.
fit_splitpop_weibull( hform, pform, data, time, event, group, start = NULL, method = "NR" )fit_splitpop_weibull( hform, pform, data, time, event, group, start = NULL, method = "NR" )
hform |
One-sided formula (e.g. '~ x1 + x2') for the H_regression equation: covariates for the log-logistic timing distribution. |
pform |
One-sided formula for the P_regression equation: covariates for the logit on the stayer ("cure") probability. |
data |
A 'data.frame' containing 'time', 'event', 'group', and every variable referenced in 'hform'/'pform'. |
time |
Character: name of the duration column ($ML_y1 in Stata). |
event |
Character: name of the 0/1 event column, 1 = failure observed, 0 = censored ($ML_y2 in Stata). |
group |
Character: name of a 0/1 column that flips which side of the cure-probability logit is used ($ML_y3 in Stata) – see the manual. If there is no such distinction in your data, pass a column of all 1s. |
start |
Optional numeric starting vector; if 'NULL', a default is constructed automatically. |
method |
Optimizer passed to [maxLik::maxLik()]; default '"NR"' (Newton-Raphson). Try '"BFGS"' if that fails to converge. |
A 'maxLik' object; use 'summary()', 'coef()', 'logLik()', 'vcov()'.
set.seed(1) n <- 150 x1 <- rnorm(n); x2 <- rbinom(n, 1, 0.5) cured <- rbinom(n, 1, plogis(0.2 + 0.5 * x2)) t_latent <- (-log(runif(n)) / exp(0.5 - 0.4 * x1))^(1 / 1.3) censor_time <- rexp(n, rate = 0.2) time <- pmax(ifelse(cured == 1, censor_time, pmin(t_latent, censor_time)), 1e-3) event <- ifelse(cured == 1, 0, as.numeric(t_latent <= censor_time)) mydata <- data.frame(time = time, event = event, x1 = x1, x2 = x2, group = 1) fit <- fit_splitpop_weibull(~x1, ~x2, mydata, "time", "event", "group", method = "BFGS") summary(fit)set.seed(1) n <- 150 x1 <- rnorm(n); x2 <- rbinom(n, 1, 0.5) cured <- rbinom(n, 1, plogis(0.2 + 0.5 * x2)) t_latent <- (-log(runif(n)) / exp(0.5 - 0.4 * x1))^(1 / 1.3) censor_time <- rexp(n, rate = 0.2) time <- pmax(ifelse(cured == 1, censor_time, pmin(t_latent, censor_time)), 1e-3) event <- ifelse(cured == 1, 0, as.numeric(t_latent <= censor_time)) mydata <- data.frame(time = time, event = event, x1 = x1, x2 = x2, group = 1) fit <- fit_splitpop_weibull(~x1, ~x2, mydata, "time", "event", "group", method = "BFGS") summary(fit)
Log-likelihood: split-population gamma model
loglik_splitpop_gamma(par, time, event, group, Xh, Xp)loglik_splitpop_gamma(par, time, event, group, Xh, Xp)
par |
Numeric parameter vector: 'H_regression' coefficients, the ancillary scalar 'theta2', then 'P_regression' coefficients. |
time, event, group
|
Numeric vectors: duration, 0/1 event indicator, 0/1 group indicator (see [fit_splitpop_loglogistic()]). |
Xh, Xp
|
Design matrices for the H_regression and P_regression parts. |
A numeric vector of per-observation log-likelihood contributions.
Xh <- cbind(1, rnorm(20)); Xp <- cbind(1, rbinom(20, 1, 0.5)) par0 <- c(1, 0, 1, 0, 0) # positive H_regression intercept: theta1 must stay > 0 loglik_splitpop_gamma(par0, time = runif(20, 0.1, 5), event = rbinom(20, 1, 0.7), group = rep(1, 20), Xh = Xh, Xp = Xp)Xh <- cbind(1, rnorm(20)); Xp <- cbind(1, rbinom(20, 1, 0.5)) par0 <- c(1, 0, 1, 0, 0) # positive H_regression intercept: theta1 must stay > 0 loglik_splitpop_gamma(par0, time = runif(20, 0.1, 5), event = rbinom(20, 1, 0.7), group = rep(1, 20), Xh = Xh, Xp = Xp)
Log-likelihood: split-population generalized gamma model
loglik_splitpop_ggamma(par, time, event, group, Xh, Xp)loglik_splitpop_ggamma(par, time, event, group, Xh, Xp)
par |
Numeric parameter vector: 'H_regression' coefficients, 'ln_sigma', 'kappa', then 'P_regression' coefficients. |
time, event, group
|
Numeric vectors: duration, 0/1 event indicator, 0/1 group indicator (see [fit_splitpop_loglogistic()]). |
Xh, Xp
|
Design matrices for the H_regression and P_regression parts. |
A numeric vector of per-observation log-likelihood contributions.
Xh <- cbind(1, rnorm(20)); Xp <- cbind(1, rbinom(20, 1, 0.5)) par0 <- c(0, 0, 0, 1, 0, 0) loglik_splitpop_ggamma(par0, time = runif(20, 0.1, 5), event = rbinom(20, 1, 0.7), group = rep(1, 20), Xh = Xh, Xp = Xp)Xh <- cbind(1, rnorm(20)); Xp <- cbind(1, rbinom(20, 1, 0.5)) par0 <- c(0, 0, 0, 1, 0, 0) loglik_splitpop_ggamma(par0, time = runif(20, 0.1, 5), event = rbinom(20, 1, 0.7), group = rep(1, 20), Xh = Xh, Xp = Xp)
Internal log-likelihood used by [fit_splitpop_loglogistic()]. Exposed so it can be evaluated by hand at candidate starting values (the R analogue of Stata's 'ml check').
loglik_splitpop_loglogistic(par, time, event, group, Xh, Xp)loglik_splitpop_loglogistic(par, time, event, group, Xh, Xp)
par |
Numeric parameter vector: 'H_regression' coefficients, the ancillary scalar 'theta2', then 'P_regression' coefficients. |
time, event, group
|
Numeric vectors: duration, 0/1 event indicator, 0/1 group indicator (see [fit_splitpop_loglogistic()]). |
Xh, Xp
|
Design matrices for the H_regression and P_regression parts. |
A numeric vector of per-observation log-likelihood contributions.
Xh <- cbind(1, rnorm(20)); Xp <- cbind(1, rbinom(20, 1, 0.5)) par0 <- c(0, 0, 1, 0, 0) loglik_splitpop_loglogistic(par0, time = runif(20, 0.1, 5), event = rbinom(20, 1, 0.7), group = rep(1, 20), Xh = Xh, Xp = Xp)Xh <- cbind(1, rnorm(20)); Xp <- cbind(1, rbinom(20, 1, 0.5)) par0 <- c(0, 0, 1, 0, 0) loglik_splitpop_loglogistic(par0, time = runif(20, 0.1, 5), event = rbinom(20, 1, 0.7), group = rep(1, 20), Xh = Xh, Xp = Xp)
Log-likelihood: split-population log-normal model
loglik_splitpop_lognormal(par, time, event, group, Xh, Xp)loglik_splitpop_lognormal(par, time, event, group, Xh, Xp)
par |
Numeric parameter vector: 'H_regression' coefficients, the ancillary scalar 'theta2', then 'P_regression' coefficients. |
time, event, group
|
Numeric vectors: duration, 0/1 event indicator, 0/1 group indicator (see [fit_splitpop_loglogistic()]). |
Xh, Xp
|
Design matrices for the H_regression and P_regression parts. |
A numeric vector of per-observation log-likelihood contributions.
Xh <- cbind(1, rnorm(20)); Xp <- cbind(1, rbinom(20, 1, 0.5)) par0 <- c(0, 0, 0, 0, 0) loglik_splitpop_lognormal(par0, time = runif(20, 0.1, 5), event = rbinom(20, 1, 0.7), group = rep(1, 20), Xh = Xh, Xp = Xp)Xh <- cbind(1, rnorm(20)); Xp <- cbind(1, rbinom(20, 1, 0.5)) par0 <- c(0, 0, 0, 0, 0) loglik_splitpop_lognormal(par0, time = runif(20, 0.1, 5), event = rbinom(20, 1, 0.7), group = rep(1, 20), Xh = Xh, Xp = Xp)
Log-likelihood: split-population Weibull model
loglik_splitpop_weibull(par, time, event, group, Xh, Xp)loglik_splitpop_weibull(par, time, event, group, Xh, Xp)
par |
Numeric parameter vector: 'H_regression' coefficients, the ancillary scalar 'theta2', then 'P_regression' coefficients. |
time, event, group
|
Numeric vectors: duration, 0/1 event indicator, 0/1 group indicator (see [fit_splitpop_loglogistic()]). |
Xh, Xp
|
Design matrices for the H_regression and P_regression parts. |
A numeric vector of per-observation log-likelihood contributions.
Xh <- cbind(1, rnorm(20)); Xp <- cbind(1, rbinom(20, 1, 0.5)) par0 <- c(0, 0, 0, 0, 0) loglik_splitpop_weibull(par0, time = runif(20, 0.1, 5), event = rbinom(20, 1, 0.7), group = rep(1, 20), Xh = Xh, Xp = Xp)Xh <- cbind(1, rnorm(20)); Xp <- cbind(1, rbinom(20, 1, 0.5)) par0 <- c(0, 0, 0, 0, 0) loglik_splitpop_weibull(par0, time = runif(20, 0.1, 5), event = rbinom(20, 1, 0.7), group = rep(1, 20), Xh = Xh, Xp = Xp)