Function for the calibration curve for the univariate joint model:
calJM <- function(object, newdata, Tstart, Thoriz, idVar = "id",
simulate = FALSE, M = 100, pl = TRUE,
include.cens = FALSE, ...) {
if (!inherits(object, "JMbayes"))
stop("Use only with 'JMbayes' objects.\n")
if (!is.data.frame(newdata) || nrow(newdata) == 0)
stop("'newdata' must be a data.frame with more than one rows.\n")
if (is.null(newdata[[idVar]]))
stop("'idVar' not in 'newdata.\n'")
id <- newdata[[idVar]]
id <- match(id, unique(id))
TermsT <- object$Terms$termsT
SurvT <- model.response(model.frame(TermsT, newdata))
is_counting <- attr(SurvT, "type") == "counting"
Time <- if (is_counting) {
ave(SurvT[, 2], id, FUN = function(x) tail(x, 1))
} else {
SurvT[, 1]
}
timeVar <- object$timeVar
newdata2 <- newdata[Time > Tstart, ]
SurvT <- model.response(model.frame(TermsT, newdata2))
if (is_counting) {
id2 <- newdata2[[idVar]]
f <- factor(id2, levels = unique(id2))
Time <- ave(SurvT[, 2], f, FUN = function(x) tail(x, 1))
delta <- ave(SurvT[, 3], f, FUN = function(x) tail(x, 1))
} else {
Time <- SurvT[, 1]
delta <- SurvT[, 2]
}
timesInd <- newdata2[[timeVar]] <= Tstart
aliveThoriz <- newdata2[Time > Thoriz & timesInd, ]
deadThoriz <- newdata2[Time <= Thoriz & delta == 1 & timesInd, ]
indCens <- Time < Thoriz & delta == 0 & timesInd
censThoriz <- newdata2[indCens, ]
nr <- length(unique(newdata2[[idVar]]))
idalive <- unique(aliveThoriz[[idVar]])
iddead <- unique(deadThoriz[[idVar]])
idcens <- unique(censThoriz[[idVar]])
Surv.aliveThoriz <- if (is_counting) {
survfitJM(object, newdata = aliveThoriz, idVar = idVar,
simulate = simulate, M = M, survTimes = Thoriz,
last.time = rep(Tstart, length(idalive)), LeftTrunc_var = all.vars(TermsT)[1L])
} else {
survfitJM(object, newdata = aliveThoriz, idVar = idVar,
simulate = simulate, M = M, survTimes = Thoriz,
last.time = rep(Tstart, length(idalive)))
}
Surv.deadThoriz <- if (is_counting) {
survfitJM(object, newdata = deadThoriz, idVar = idVar,
simulate = simulate, survTimes = Thoriz,
last.time = rep(Tstart, length(iddead)), LeftTrunc_var = all.vars(TermsT)[1L])
} else {
survfitJM(object, newdata = deadThoriz, idVar = idVar,
simulate = simulate, survTimes = Thoriz, last.time = rep(Tstart, length(iddead)))
}
Surv.aliveThoriz <- sapply(Surv.aliveThoriz$summaries, "[", 2)
Surv.deadThoriz <- sapply(Surv.deadThoriz$summaries, "[", 2)
if (nrow(censThoriz)) {
Surv.censThoriz <- if (is_counting) {
survfitJM(object, newdata = censThoriz, idVar = idVar,
simulate = simulate, M = M, survTimes = Thoriz,
last.time = rep(Tstart, length(idcens)), LeftTrunc_var = all.vars(TermsT)[1L])
} else {
survfitJM(object, newdata = censThoriz, idVar = idVar,
simulate = simulate, M = M, survTimes = Thoriz,
last.time = rep(Tstart, length(idcens)))
}
tt <- Time[indCens]
weights <- if (is_counting) {
survfitJM(object, newdata = censThoriz, idVar = idVar,
simulate = simulate, M = M, survTimes = Thoriz,
last.time = tt[!duplicated(censThoriz[[idVar]])],
LeftTrunc_var = all.vars(TermsT)[1L])
}
else {
survfitJM(object, newdata = censThoriz, idVar = idVar,
simulate = simulate, M = M, survTimes = Thoriz,
last.time = tt[!duplicated(censThoriz[[idVar]])])
}
Surv.censThoriz <- sapply(Surv.censThoriz$summaries,
"[", 2)
weights <- sapply(weights$summaries, "[", 2)
} else {
Surv.censThoriz <- weights <- NA
}
# until here I used the prederrJM function from the JMbayes package
newdata.id <- newdata[tapply(row.names(newdata), newdata$id, tail, 1), ]
statusCens <- rbinom(n = length(weights), size = 1, prob = weights)
matCens <- data.frame(years = newdata.id$years[newdata.id$id %in% names(Surv.censThoriz)],
pred_prob = 1 - Surv.censThoriz, status = 1 - statusCens)
# New part for calibration plot
matDead <- data.frame(years = newdata.id$years[newdata.id$id %in% names(Surv.deadThoriz)],
pred_prob = 1 - Surv.deadThoriz, status = 1)
matAlive <- data.frame(years = newdata.id$years[newdata.id$id %in% names(Surv.aliveThoriz)],
pred_prob = 1 - Surv.aliveThoriz, status = 0)
if (include.cens == TRUE) {
mat <- rbind(matDead, matAlive, matCens)
} else {
mat <- rbind(matDead, matAlive)
}
# Cox model calibration curve
# The complementary log-log transformation for the predicted probabilities
# rather than the probabilities themselves.
mat$cox.cll <- log(-log(1 - mat$pred_prob))
# Use a conventional Cox proportional hazards model with restricted cubic
# splines to model the relationship between
# log(-log(1-pred)) and the log-hazard of the outcome
calibrate.cox <- coxph(Surv(years, status) ~ rcs(cox.cll, 3), x = T,
data = mat)
# Based on the fitted model, an estimated probability of the occurrence of the
# outcome prior to time t_0 can be estimated for each value of pred
predict.grid.cox <- seq(0.01,
0.99, length = 100)
predict.grid.cox.cll <- log(-log(1 - predict.grid.cox))
predict.grid.cox.df <- data.frame(predict.grid.cox)
predict.grid.cox.cll.df <- data.frame(predict.grid.cox.cll)
names(predict.grid.cox.df) <- "cox"
names(predict.grid.cox.cll.df) <- "cox.cll"
predict.calibrate.cox <- 1 - predictSurvProb(calibrate.cox,
newdata = predict.grid.cox.cll.df,
times = Thoriz)
# Save results
out <- list(x = predict.grid.cox, y = predict.calibrate.cox)
if (pl == TRUE) {
plot(predict.grid.cox, predict.calibrate.cox, type="l", lty=1, col="red",
xlim=c(0,1),ylim=c(0,1),
xlab = "Predicted probability",
ylab = "Observed probability")
abline(0, 1, lwd = 6, col = gray(0.85))
}
out
}Function for the calibration curve for the multivariate joint model:
calJM.mv <- function(object, newdata, Tstart, Thoriz, idVar = "id",
simulate = FALSE, M = 100, pl = TRUE,
include.cens = FALSE, ...) {
if (!inherits(object, "mvJMbayes"))
stop("Use only with 'mvJMbayes' objects.\n")
if (!is.data.frame(newdata) || nrow(newdata) == 0)
stop("'newdata' must be a data.frame with more than one rows.\n")
if (is.null(newdata[[idVar]]))
stop("'idVar' not in 'newdata.\n'")
id <- newdata[[idVar]]
id <- match(id, unique(id))
TermsT <- object$model_info$coxph_components$Terms
environment(TermsT) <- parent.frame()#.GlobalEnv
SurvT <- model.response(model.frame(TermsT, newdata))
is_counting <- attr(SurvT, "type") == "counting"
is_interval <- attr(SurvT, "type") == "interval"
Time <- if (is_counting) {
ave(SurvT[, 2], id, FUN = function (x) tail(x, 1))
} else if (is_interval) {
Time1 <- SurvT[, "time1"]
Time2 <- SurvT[, "time2"]
Time <- Time1
Time[Time2 != 1] <- Time2[Time2 != 1]
Time
} else {
SurvT[, 1]
}
timeVar <- object$model_info$timeVar
newdata2 <- newdata[Time > Tstart, ]
id2 <- newdata2[[idVar]]
SurvT <- model.response(model.frame(TermsT, newdata2))
if (is_counting) {
f <- factor(id2, levels = unique(id2))
Time <- ave(SurvT[, 2], f, FUN = function (x) tail(x, 1))
event <- ave(SurvT[, 3], f, FUN = function (x) tail(x, 1))
} else if (is_interval) {
Time1 <- SurvT[, "time1"]
Time2 <- SurvT[, "time2"]
Time <- Time1
Time[Time2 != 1] <- Time2[Time2 != 1]
event <- SurvT[, "status"]
} else {
Time <- SurvT[, 1]
event <- SurvT[, 2]
}
timesInd <- newdata2[[timeVar]] <= Tstart
aliveThoriz <- newdata2[Time > Thoriz & timesInd, ]
deadThoriz <- newdata2[Time <= Thoriz & (event == 1 | event == 3) & timesInd, ]
indCens <- Time < Thoriz & (event == 0 | event == 2) & timesInd
censThoriz <- newdata2[indCens, ]
nr <- length(unique(newdata2[[idVar]]))
idalive <- unique(aliveThoriz[[idVar]])
iddead <- unique(deadThoriz[[idVar]])
idcens <- unique(censThoriz[[idVar]])
Surv.aliveThoriz <- if (is_counting) {
survfitJM(object, newdata = aliveThoriz, idVar = idVar, M = M,
survTimes = Thoriz, last.time = rep(Tstart, length(idalive)),
LeftTrunc_var = all.vars(TermsT)[1L])
} else {
survfitJM(object, newdata = aliveThoriz, idVar = idVar, M = M,
survTimes = Thoriz, last.time = rep(Tstart, length(idalive)))
}
Surv.deadThoriz <- if (is_counting) {
survfitJM(object, newdata = deadThoriz, idVar = idVar,
survTimes = Thoriz, last.time = rep(Tstart, length(iddead)),
LeftTrunc_var = all.vars(TermsT)[1L])
} else {
survfitJM(object, newdata = deadThoriz, idVar = idVar,
survTimes = Thoriz, last.time = rep(Tstart, length(iddead)))
}
Surv.aliveThoriz <- sapply(Surv.aliveThoriz$summaries, "[", 2)
Surv.deadThoriz <- sapply(Surv.deadThoriz$summaries, "[", 2)
if (nrow(censThoriz)) {
Surv.censThoriz <- if (is_counting) {
survfitJM(object, newdata = censThoriz, idVar = idVar, M = M,
survTimes = Thoriz, last.time = rep(Tstart, length(idcens)),
LeftTrunc_var = all.vars(TermsT)[1L])
} else {
survfitJM(object, newdata = censThoriz, idVar = idVar, M = M,
survTimes = Thoriz, last.time = rep(Tstart, length(idcens)))
}
tt <- Time[indCens]
weights <- if (is_counting) {
survfitJM(object, newdata = censThoriz, idVar = idVar, M = M,
survTimes = Thoriz, last.time = tt[!duplicated(censThoriz[[idVar]])],
LeftTrunc_var = all.vars(TermsT)[1L])
} else {
survfitJM(object, newdata = censThoriz, idVar = idVar, M = M,
survTimes = Thoriz, last.time = tt[!duplicated(censThoriz[[idVar]])])
}
Surv.censThoriz <- sapply(Surv.censThoriz$summaries, "[", 2)
weights <- sapply(weights$summaries, "[", 2)
} else {
Surv.censThoriz <- weights <- NA
}
# until here I used the prederrJM.mv function from the JMbayes package
newdata.id <- newdata[tapply(row.names(newdata), newdata$id, tail, 1), ]
statusCens <- rbinom(n = length(weights), size = 1, prob = weights)
matCens <- data.frame(years = newdata.id$years[newdata.id$id %in% names(Surv.censThoriz)],
pred_prob = 1 - Surv.censThoriz, status = 1 - statusCens)
# New part for calibration plot
matDead <- data.frame(years = newdata.id$years[newdata.id$id %in% names(Surv.deadThoriz)],
pred_prob = 1 - Surv.deadThoriz, status = 1)
matAlive <- data.frame(years = newdata.id$years[newdata.id$id %in% names(Surv.aliveThoriz)],
pred_prob = 1 - Surv.aliveThoriz, status = 0)
if (include.cens == TRUE) {
mat <- rbind(matDead, matAlive, matCens)
} else {
mat <- rbind(matDead, matAlive)
}
# Cox model calibration curve
# The complementary log-log transformation for the predicted probabilities
# rather than the probabilities themselves.
mat$cox.cll <- log(-log(1 - mat$pred_prob))
# Use a conventional Cox proportional hazards model with restricted cubic
# splines to model the relationship between
# log(-log(1-pred)) and the log-hazard of the outcome
calibrate.cox <- coxph(Surv(years, status) ~ rcs(cox.cll, 3), x = T,
data = mat)
# Based on the fitted model, an estimated probability of the occurrence of the
# outcome prior to time t_0 can be estimated for each value of pred
predict.grid.cox <- seq(0.01,
0.99, length = 100)
predict.grid.cox.cll <- log(-log(1 - predict.grid.cox))
predict.grid.cox.df <- data.frame(predict.grid.cox)
predict.grid.cox.cll.df <- data.frame(predict.grid.cox.cll)
names(predict.grid.cox.df) <- "cox"
names(predict.grid.cox.cll.df) <- "cox.cll"
predict.calibrate.cox <- 1 - predictSurvProb(calibrate.cox,
newdata = predict.grid.cox.cll.df,
times = Thoriz)
# Save results
out <- list(x = predict.grid.cox, y = predict.calibrate.cox)
if (pl == TRUE) {
plot(predict.grid.cox, predict.calibrate.cox, type="l", lty=1, col="red",
xlim=c(0,1),ylim=c(0,1),
xlab = "Predicted probability",
ylab = "Observed probability")
abline(0, 1, lwd = 6, col = gray(0.85))
}
out
}© Eleni-Rosalina Andrinopoulou