Nothing
intcure <- function(formula = formula(data), cureform = NULL,
data = sys.parent(), subset, na.action, bt = NULL, gm = NULL,
basepara = NULL, sigma = c(0, 0, 0),
optimcfg = list(ndeps = 0.001, maxit = 1000, reltol = 0.00001,
method = "Nelder-Mead", hessian = FALSE),
intcfg = list(eps = 0.0001, lower = c(-5, -5), upper = c(5, 5)),
basedist = c("exponential", "weibull", "piecewise"),
npiece = 3, piececut = NULL, piececuttype = c("quantile", "even"),
model = FALSE, y = TRUE, x = FALSE, z = FALSE, funval = FALSE,
debug = c("intcure", "integration", "optim"))
{
call <- match.call()
formal.args <- formals(sys.function(sys.parent()))
m <- match.call(expand.dots = FALSE)
temp <- c("", "formula", "data", "subset", "na.action")
m1 <- m[match(temp, names(m), nomatch = 0)]
temp <- c("", "cureform", "data", "subset", "na.action")
m2 <- m[match(temp, names(m), nomatch = 0)]
if (missing(debug)) debug <- NULL
else debug <- match.arg(debug, several.ok = TRUE)
special <- "cluster"
Terms1 <- if(missing(data)) terms(formula, special)
else terms(formula, special, data = data)
m1$formula <- Terms1
m1[[1]] <- as.name("model.frame")
if(m1[[2]][[1]] == as.name("Surv")) {
m1[[2]][[2]] <- m1[[2]]
m1[[2]][[1]] <- as.name("~")
m1[[2]][[3]] <- 1
}
m1 <- eval(m1, sys.parent())
rawY <- model.extract(m1, "response")
if(!inherits(rawY, "Surv"))
stop("Response must be a survival object")
sorted <- order(rawY[, 1], 1-rawY[,2])
Y <- rawY[sorted, ]
datalen <- nrow(Y)
offsetx <- attr(Terms1, "offset")
tt <- length(offsetx)
offsetx <- if(tt == 0) rep(0, datalen) else if(tt == 1)
m1[[offsetx]]
else {
ff <- m1[[offsetx[1]]]
for(i in 2:tt)
ff <- ff + m1[[offsetx[i]]]
ff
}
offsetx <- offsetx[sorted]
clusterx <- attr(Terms1, "specials")$cluster
dropx <- NULL
nclusterx <- 0
if(length(clusterx)) {
tempc <- untangle.specials(Terms1, "cluster", 1:10)
ord <- attr(Terms1, "order")[tempc$terms]
if(any(ord > 1))
stop("Cluster can not be used in an interaction")
clusterx <- strata(m1[, tempc$vars], shortlabel = TRUE)
dropx <- tempc$terms
clusterx <- as.numeric(clusterx)[sorted]
nclusterx <- length(unique(clusterx))
}
else {
clusterx <- 1:datalen
nclusterx <- datalen
}
if(length(dropx))
X <- model.matrix(Terms1[ - dropx], m1)[, -1, drop = FALSE]
else X <- model.matrix(Terms1, m1)[, -1, drop = FALSE]
X <- X[sorted, , drop = FALSE]
nbt <- ncol(X)
if (!length(X)) X <- NULL
if (is.null(nbt)) nbt <- 0
cofnames <- dimnames(X)[[2]]
ngm <- 0
Z <- NULL
clusterz <- NULL
nclusterz <- 0
offsetz <- NULL
if (!missing(cureform)) {
Terms2 <- if(missing(data)) terms(cureform, special)
else terms(cureform, special, data = data)
m2$formula <- Terms2
m2$cureform <- NULL
m2[[1]] <- as.name("model.frame")
curedata <- eval(m2, sys.parent())
Terms2 <- attr(curedata, "terms")
offsetz <- attr(Terms2, "offset")
tt <- length(offsetz)
offsetz <- if(tt == 0) rep(0, datalen) else if(tt == 1)
curedata[[offsetz]]
else {
ff <- curedata[[offsetz[1]]]
for(i in 2:tt)
ff <- ff + curedata[[offsetz[i]]]
ff
}
offsetz <- offsetz[sorted]
clusterz <- attr(Terms2, "specials")$cluster
dropz <- NULL
if(length(clusterz)) {
tempc <- untangle.specials(Terms2, "cluster", 1:10)
ord <- attr(Terms2, "order")[tempc$terms]
if(any(ord > 1))
stop("Cluster can not be used in an interaction")
clusterz <- strata(curedata[, tempc$vars], shortlabel = TRUE)
dropz <- tempc$terms
clusterz <- as.numeric(clusterz)[sorted]
}
else {
clusterz <- 1:datalen
nclusterz <- datalen
}
if(length(dropz))
Z <- model.matrix(Terms2[ - dropz], curedata)[, , drop = FALSE]
else Z <- model.matrix(Terms2, curedata)[, , drop = FALSE]
Z <- Z[sorted, , drop = FALSE]
ngm <- ncol(Z)
cofnames <- c(cofnames,dimnames(Z)[[2]])
nclusterz <- length(unique(clusterz))
}
basedist <- match.arg(basedist)
piececuttype <- match.arg(piececuttype)
names(sigma) = c("lsigmau", "lsigmav", "ftrho")
if (any(is.na(sigma[1:2])) && !is.na(sigma[3]))
stop("If the covariance is present, the two random effects must both be present.")
intcfg.def = eval(formal.args[["intcfg"]], envir = sys.frame(sys.parent()))
missed = names(intcfg.def)[!(names(intcfg.def) %in% names(intcfg))]
intcfg[missed] = intcfg.def[missed]
if (!is.null(debug) && any(debug=="integration")) intcfg$trace = TRUE
else intcfg$trace = FALSE
optimcfg.def = eval(formal.args[["optimcfg"]], envir = sys.frame(sys.parent()))
missed = names(optimcfg.def)[!(names(optimcfg.def) %in% names(optimcfg))]
optimcfg[missed] = optimcfg.def[missed]
if (!is.null(debug) && any(debug=="optim")) optimcfg$trace = 2
else optimcfg$trace = 0
ptime = NULL
if (basedist=="piecewise") {
pdist = piecetime(Y, piececuttype, piececut, npiece)
nbasepara = pdist$nbasepara
ptime = pdist$ptime
}
else if (basedist=="exponential") nbasepara = 1
else nbasepara = 2
if(model) {
list(Y = Y, Z = Z, X = X, ngm = ngm, nbt = nbt, nbasepara = nbasepara,
ptime = ptime)
}
else {
para = c(initpara(formula, cureform, basedist, bt, gm, basepara, nbt,
ngm, nbasepara, data, rawY), list(sigma = sigma))
data2 <- list(y = Y, x = X, z = Z, id = clusterx, ncluster = nclusterx,
offsetx = offsetx, offsetz = offsetz)
if (!is.null(debug) && any(debug=="intcure")) {
fcat("Initial values:\n")
fcat("sigma =", para$sigma, "\n")
fcat("bt =", para$bt, "\n")
fcat("gm =", para$gm, "\n")
fcat("basepara =", para$basepara, "\n")
}
fit <- intcure.fit(para, data2, intcfg, optimcfg, basedist, funval, ptime)
if (is.null(names(fit$bt))) names(fit$bt) <- dimnames(X)[[2]]
if (is.null(names(fit$gm))) names(fit$gm) <- dimnames(Z)[[2]]
fit$call <- call
fit$basedist <- basedist
fit$ptime <- ptime
fit$n <- datalen
class(fit) <- "intcure"
fit
}
}
parv2l <- function(vpara, lpara)
{
nbt = length(lpara$bt)
ngm = length(lpara$gm)
nbasepara = length(lpara$basepara)
np = length(vpara)
if (nbt > 0) bt = vpara[1:nbt]
else bt = NULL
gm = vpara[nbt+(1:ngm)]
basepara = vpara[nbt + ngm + (1:nbasepara)]
sigma = lpara$sigma
if (np > nbt + ngm + nbasepara)
sigma[!is.na(sigma)] = vpara[(nbt + ngm + nbasepara + 1):np]
list(bt = bt, gm = gm, basepara = basepara, sigma = sigma)
}
piecetime = function(Y, piececuttype, piececut, npiece)
{
maxtime <- max(Y[Y[, 2]==1, 1])
mintime <- min(Y[Y[, 2]==1, 1])
if (is.null(piececut)) {
nbasepara = npiece
if (piececuttype == "even")
ptime <- mintime + (1:(npiece-1))*(maxtime - mintime)/npiece
else ptime <- quantile(Y[Y[, 2]==1, 1],
probs = seq(0, 1, 1/npiece))[-c(1, npiece+1)]
}
else {
if (max(piececut) >= maxtime || min(piececut) <= mintime)
stop("piececut values are beyond survival time range.")
nbasepara = length(piececut) + 1
ptime <- piececut
}
list(ptime = ptime, nbasepara = nbasepara)
}
baseline <- function(basepara, time, basedist, ptime)
{
nbasepara = length(basepara)
if (basedist == "exponential") {
rate <- exp(basepara)
h <- rep(rate, length(time))
S <- exp(-rate*time)
list(h = h, S = S)
}
else if (basedist == "weibull") {
p <- exp(basepara[1])
rate <- exp(basepara[2])
h <- rate*p*time^(p-1.0)
S <- exp(-rate*time^p)
list(h = h, S = S)
}
else if (basedist == "piecewise") {
hazard <- exp(basepara)
num <- as.numeric(cut(time, breaks = c(0, ptime, Inf), right = FALSE,
include.lowest = TRUE))
h <- hazard[num]
cumhaz <- c(0, cumsum(hazard[-nbasepara]*diff(c(0, ptime))))
chaz <- h*(time - c(0, ptime[-nbasepara])[num])
S <- exp(-(cumhaz[num] + chaz))
list(h = h, S = S)
}
else stop("baseline distribution is not defined.")
}
basesurvfun = function(time, intcureobj)
{
basedist = intcureobj$basedist
basepara = intcureobj$basepara
nbasepara = length(basepara)
if (basedist == "exponential") {
rate <- exp(basepara)
exp(-rate*time)
}
else if (basedist == "weibull") {
p <- exp(basepara[1])
rate <- exp(basepara[2])
exp(-rate*time^p)
}
else if (basedist == "piecewise") {
ptime = intcureobj$ptime
hazard <- exp(basepara)
num <- as.numeric(cut(time, breaks = c(0, ptime), right = FALSE,
include.lowest = TRUE))
cumhaz <- c(0, cumsum(hazard*diff(c(0,
ptime)))[-nbasepara])
chaz <- hazard[num]*(time - c(0, ptime[-nbasepara])[num])
exp(-(cumhaz[num] + chaz))
}
else stop("baseline distribution is not defined.")
}
clusterlik <- function(bx, gz, basesurv, basehaz, cens)
{
uncure <- exp(gz)/(1+exp(gz))
surv <- basesurv^(exp(bx))
val <- sum(ifelse(cens==1, gz - log(1.0 + exp(gz))+log(basehaz) + bx +
exp(bx)*log(basesurv),
-log(1.0 + exp(gz)) + log(1.0 + exp(gz)*surv)))
val
}
clusterliku <- function(u, bx, gz, basesurv, basehaz, cens, sigma)
{
sapply(u, function(x, bx, gz, basesurv, basehaz, cens, sigma){
uncure <- exp(gz)/(1+exp(gz))
surv <- basesurv^(exp(bx+x))
exp(sum(ifelse(cens==1, gz - log(1+exp(gz)) +
log(basehaz) + bx + x + exp(bx+x)*log(basesurv),
log(1-uncure+uncure*surv))) + dnorm(x, sd = exp(sigma), log = TRUE))
}, bx, gz, basesurv, basehaz, cens, sigma)
}
clusterlikv <- function(v, bx, gz, basesurv, basehaz, cens, sigma)
{
sapply(v, function(x, bx, gz, basesurv, basehaz, cens, sigma){
uncure <- exp(gz + x)/(1+exp(gz + x))
surv <- basesurv^(exp(bx))
exp(sum(ifelse(cens==1, gz + x - log(1+exp(gz + x)) +
log(basehaz) + bx + exp(bx)*log(basesurv),
log(1-uncure+uncure*surv))) + dnorm(x, sd = exp(sigma), log = TRUE))
}, bx, gz, basesurv, basehaz, cens, sigma)
}
clusterlikuv <- function(uv, bx, gz, basesurv, basehaz, cens, sigma)
{
uncure <- exp(gz + uv[2])/(1+exp(gz + uv[2]))
surv <- basesurv^(exp(bx+uv[1]))
val <- sum(ifelse(cens==1, gz + uv[2] - log(1+exp(gz + uv[2])) +
log(basehaz) + bx+uv[1] + exp(bx+uv[1])*log(basesurv),
log(1-uncure+uncure*surv)))
tmp <- dnorm(uv[1], sd = exp(sigma[1]), log = TRUE) +
dnorm(uv[2], sd = exp(sigma[2]), log = TRUE)
exp(val + tmp)
}
clusterlikuv2 <- function(uv, bx, gz, basesurv, basehaz, cens, sigma)
{
uncure <- exp(gz + uv[2])/(1+exp(gz + uv[2]))
surv <- basesurv^(exp(bx+uv[1]))
rho = (exp(2.0*sigma[3])-1.0)/(exp(2.0*sigma[3])+1.0)
cvar <- rho*exp(sigma[1]+sigma[2])
val <- sum(ifelse(cens==1, gz + uv[2] - log(1+exp(gz + uv[2])) +
log(basehaz) + bx+uv[1] + exp(bx+uv[1])*log(basesurv),
log(1-uncure+uncure*surv)))
tmp <- dmvnorm(uv, sigma = matrix(c(exp(2.0*sigma[1]), cvar, cvar,
exp(2.0*sigma[2])), nrow = 2), log = TRUE)
exp(val + tmp)
}
loglikfun <- function(vpara, data, basedist, intcfg, lpara, failinfo)
{
pars <- parv2l(vpara, lpara)
time <- data$y[, 1]
if (is.null(data$x) || length(pars$bt)==0) bx = rep(0, length(time))
else bx <- data$x %*% pars$bt
gz <- data$z %*% pars$gm
sigma <- pars$sigma
hS <- baseline(pars$basepara, time, basedist, failinfo)
if (intcfg$trace) fcat("Parameter =", vpara, "\n")
val <- 0.0
for (i in 1:data$ncluster) {
member <- data$id==i
cens <- data$y[member, 2]
basesurv <- hS$S[member]
basehaz <- hS$h[member]
cbx <- bx[member]
cgz <- gz[member]
if (all(is.na(lpara$sigma))) {
intv = list()
intv$value = clusterlik(cbx, cgz, basesurv, basehaz, cens)
intv$error <- intv$code <- NULL
}
else if (all(is.na(lpara$sigma)[2:3])) {
intcfg$fun = "integrate"
intv <- intfun(clusterliku, intcfg, bx = cbx, gz = cgz,
basesurv = basesurv, basehaz = basehaz, cens = cens,
sigma = sigma[1])
}
else if (all(is.na(lpara$sigma)[c(1, 3)])) {
intcfg$fun = "integrate"
intv <- intfun(clusterlikv, intcfg, bx = cbx, gz = cgz,
basesurv = basesurv, basehaz = basehaz, cens = cens,
sigma = sigma[2])
}
else if (is.na(lpara$sigma[3])) {
if(is.null(intcfg$fun)) intcfg$fun = "adaptIntegrate"
intv <- intfun(clusterlikuv, intcfg, bx = cbx, gz = cgz,
basesurv = basesurv, basehaz = basehaz, cens = cens,
sigma = sigma[1:2])
}
else {
if(is.null(intcfg$fun)) intcfg$fun = "adaptIntegrate"
intv <- intfun(clusterlikuv2, intcfg, bx = cbx, gz = cgz,
basesurv = basesurv, basehaz = basehaz, cens = cens,
sigma = sigma)
}
if (intcfg$trace)
fcat("Cluster", i, "likelihood =", intv$value, "with code =",
intv$code, "and error =", intv$error, "\n")
val <- val + intv$value
}
if (intcfg$trace) fcat("total likelihood =", val, "\n")
val
}
intfun = function(fn, intcfg, bx, gz, basesurv, basehaz, cens, sigma)
{
if (intcfg$fun == "integrate") {
intv = integrate(fn, lower = intcfg$lower[1], upper = intcfg$upper[1],
bx = bx, gz = gz, basesurv = basesurv, basehaz = basehaz,
cens = cens, sigma = sigma, rel.tol = intcfg$eps)
list(value = log(intv$value), error = intv$abs.error, code = intv$message)
}
else if (intcfg$fun == "adaptIntegrate") {
if (!exists("adaptIntegrate"))
stop("adaptIntegrate cannot be found. It requires cubature package\n")
intv <- adaptIntegrate(fn, lowerLimit = intcfg$lower,
upperLimit = intcfg$upper, tol = intcfg$eps, bx = bx, gz = gz,
basesurv = basesurv, basehaz = basehaz, cens = cens, sigma = sigma)
list(value = log(intv$integral), error = intv$error, code = intv$returnCode)
}
else stop("Unknown integration function in intcfg$fun\n")
}
intcure.fit <- function(lpara, data, intcfg, optimcfg, basedist, funval, failinfo)
{
vpara <- c(lpara$bt, lpara$gm, lpara$basepara, lpara$sigma[!is.na(lpara$sigma)])
if (funval) {
val <- lpara
val$loglik <- loglikfun(vpara, data = data, basedist = basedist,
intcfg = intcfg, lpara = lpara, failinfo = failinfo)
val
}
else {
method = optimcfg$method
hessian = optimcfg$hessian
optimcfg$method <- optimcfg$hessian <- NULL
if (length(optimcfg$ndeps)==1)
optimcfg$ndeps = rep(optimcfg$ndeps, length(vpara))
else {
if (length(optimcfg$ndeps)!=length(vpara))
stop("ndeps is of the wrong length")
}
fit <- optim(vpara, loglikfun, method = method,
control = c(list(fnscale = -1), optimcfg), data = data,
basedist = basedist, intcfg = intcfg, lpara = lpara,
failinfo = failinfo, hessian = hessian)
val <- parv2l(fit$par, lpara)
if (all(is.na(lpara$sigma))) val$sigma = NULL
else val$sigma = val$sigma[(1:3)[!is.na(lpara$sigma)]]
val$loglik <- fit$value
if (hessian) val$hessian <- fit$hessian
val
}
}
print.intcure <- function(x, ...)
{
if(!is.null(cl <- x$call)) {
cat("Call:\n")
dput(cl)
}
cat("\nMixture cure model with", x$basedist, "baseline\n")
cat("\nFixed effects in survival model:\n")
print(x$bt, ...)
cat("\nBaseline parameters in survival model:\n")
print(x$basepara, ...)
cat("\nFixed effects in logistic model:\n")
print(x$gm, ...)
cat("\nRandom effects:\n")
print(x$sigma, ...)
cat("\nMaximum log-likelihood:", x$loglik, "\n")
invisible(x)
}
summary.intcure = function(object, ...)
{
if (is.null(object$hessian)) stop("Hessian matrix is needed")
varpar = solve(-object$hessian)
stderr = parv2l(sqrt(diag(varpar)), object[c("bt", "gm", "basepara", "sigma")])
object$statistics = list(
bt = cbind(Estimate = object$bt, Stderr = stderr$bt,
zvalue = object$bt/stderr$bt,
pvalue = (1-pnorm(abs(object$bt/stderr$bt)))*2),
gm = cbind(Estimate = object$gm, Stderr = stderr$gm,
zvalue = object$gm/stderr$gm,
pvalue = (1-pnorm(abs(object$gm/stderr$gm)))*2),
basepara = cbind(Estimate = object$basepara, Stderr = stderr$basepara,
zvalue = object$basepara/stderr$basepara,
pvalue = (1-pnorm(abs(object$basepara/stderr$basepara)))*2),
sigma = cbind(Estimate = object$sigma, Stderr = stderr$sigma))
class(object) = "summary.intcure"
object
}
print.summary.intcure = function(x, ...)
{
if(!is.null(cl <- x$call)) {
cat("Call:\n")
dput(cl)
}
cat("\nMixture cure model with", x$basedist, "baseline\n")
cat("\nFixed effects in survival model:\n")
print(x$statistics$bt, ...)
cat("\nBaseline parameters in survival model:\n")
print(x$statistics$basepara, ...)
cat("\nFixed effects in logistic model:\n")
print(x$statistics$gm, ...)
cat("\nRandom effects:\n")
print(x$statistics$sigma, ...)
cat("\nMaximum log-likelihood:", x$loglik, "\n")
invisible(x)
}
coef.intcure <- function(object, ...)
{
c(object$bt, object$gm, object$basepara, object$sigma)
}
loglik.intcure <- function(x)x$loglik
initpara = function(formula, cureform, basedist, bt, gm, basepara, nbt, ngm,
nbasepara, data, Y)
{
phcall <- gsub(" *[+] *cluster[(][^)]*[)] *[+]*", "",
paste(deparse(formula, width.cutoff = 500), collapse = ""))
logitcall <- gsub(" *[+] *cluster[(][^)]*[)] *[+]*", "",
paste(deparse(cureform, width.cutoff = 500), collapse = ""))
lfit = survreg(eval(parse(text = phcall)), data = data,
dist = "exponential", subset = Y[, 1] <= max(Y[Y[, 2]==1, 1]))
if (is.null(bt)) bt <- coef(lfit)[-1]
else names(bt) = names(coef(lfit)[-1])
maxtime = max(Y[Y[, 2]==1, 1])
mintime = min(Y[Y[, 2]==1, 1])
uncure = ifelse(Y[, 2] == 1 | Y[, 1] < mintime, 1,
ifelse(Y[, 1] > maxtime, 0,
1-punif(Y[, 1], min = mintime, max = maxtime)))
ifit = glm(eval(parse(text = paste("uncure", logitcall))),
data = cbind(uncure, data), family=binomial(link=logit))
if (is.null(gm)) gm <- coef(ifit)
else names(gm) = names(coef(ifit))
if (length(bt)!=nbt || length(gm)!=ngm)
stop("The length of initial values doesn't match needs.")
if (is.null(basepara)) {
if (basedist=="exponential")
basepara <- c(lograte = unname(-coef(lfit)[1]))
else if (basedist=="weibull")
basepara <- c(logshape = 0, lograte = unname(-coef(lfit)[1]))
else if (basedist=="piecewise")
basepara <- structure(rep(-coef(lfit)[1], nbasepara),
names = paste("loghazard", 1:nbasepara, sep = ""))
else stop("basedist is not valid in baseinit.")
}
else {
if (basedist=="exponential") names(basepara) <- "lograte"
else if (basedist=="weibull") names(basepara) <- c("logshape",
"lograte")
else if (basedist=="piecewise")
names(basepara) <- paste("loghazard", 1:nbasepara, sep = "")
else stop("basedist is not valid in baseinit.")
}
if (length(bt) != nbt)
stop("The number of regression parameters in latency is not correct.")
if (length(gm) != ngm)
stop("The number of regression parameters in incidence is not correct.")
if (length(basepara) != nbasepara)
stop("The number of baseline parameters is not correct.")
list(bt = bt, gm = gm, basepara = basepara)
}
fcat = function (...)
{
cat(...)
flush.console()
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.