Nothing
##' @title Prediction in Joint Models
##' @name survfitJMMLSM
##' @aliases survfitJMMLSM
##' @description This function computes the conditional probability of
##' surviving later times than the last observed time for which a longitudinal
##' measurement was available.
##' @param object an object inheriting from class \code{JMMLSM}.
##' @param seed a random seed number to proceed non-parametric bootstrap. Default is 100.
##' @param ynewdata a data frame that contains the longitudinal and covariate information for the subjects
##' for which prediction of survival probabilities is required.
##' @param cnewdata a data frame that contains the survival and covariate information for the subjects
##' for which prediction of survival probabilities is required.
##' @param u a numeric vector of times for which prediction survival probabilities are to be computed.
##' @param Last.time a numeric vector or character string. This specifies the known time at which each of
##' the subjects in cnewdata was known to be alive. If NULL, then this is automatically taken as the
##' survival time of each subject. If a numeric vector, then it is assumed to be greater than or equals to the
##' last available longitudinal time point for each subject. If a character string, then it should be
##' a variable in cnewdata.
##' @param obs.time a character string of specifying a longitudinal time variable in ynewdata.
##' @param method a character string specifying the type of probability approximation; if \code{Laplace}, then a first order estimator is computed.
##' If \code{GH}, then the standard Gauss-Hermite quadrature is used instead.
##' @param quadpoint number of quadrature points used for estimating conditional probabilities
##' when \code{method = "GH"}. Default is NULL. If \code{method = "GH"}, then 15 is used.
##' @param ... further arguments passed to or from other methods.
##' @return a list of matrices with conditional probabilities for subjects.
##' @author Shanpeng Li \email{lishanpeng0913@ucla.edu}
##' @seealso \code{\link{JMMLSM}}
##' @export
##'
survfitJMMLSM <- function(object, seed = 100, ynewdata = NULL, cnewdata = NULL,
u = NULL, Last.time = NULL, obs.time = NULL, method = c("Laplace", "GH"), quadpoint = NULL, ...) {
if (!inherits(object, "JMMLSM"))
stop("Use only with 'JMMLSM' objects.\n")
if (is.null(ynewdata))
stop("New longitudinal data for dynamic prediction is needed.")
if (is.null(cnewdata))
stop("New longitudinal data for dynamic prediction is needed.")
if (is.null(u))
stop("Please specify the future time for dynamic prediction.")
if (!method %in% c("Laplace", "GH"))
stop("Please specify a method for probability approximation: Laplace or GH.")
if (!is.vector(u))
stop("u must be vector typed.")
if (is.null(quadpoint)) {
quadpoint <- object$quadpoint
}
if (is.null(obs.time)) {
stop("Please specify a vector that represents the time variable from ydatanew.")
} else {
if (!obs.time %in% colnames(ynewdata)) {
stop(paste0(obs.time, " is not found in ynewdata."))
}
}
bvar <- all.vars(object$random)
ID <- bvar[length(bvar)]
if (!(ID %in% colnames(ynewdata)))
stop(paste("The ID variable", ID, "is not found in ynewdata."))
if (!(ID %in% colnames(cnewdata)))
stop(paste("The ID variable", ID, "is not found in cnewdata."))
ynewdata <- ynewdata[, colnames(object$ydata)]
cnewdata <- cnewdata[, colnames(object$cdata)]
ydata2 <- rbind(object$ydata, ynewdata)
cdata2 <- rbind(object$cdata, cnewdata)
variance.formula <- as.formula(paste("", object$LongitudinalSubmodelvariance[3], sep = "~"))
getdum <- getdummy(long.formula = object$LongitudinalSubmodelmean,
surv.formula = object$SurvivalSubmodel,
variance.formula = variance.formula,
random = object$random, ydata = ydata2, cdata = cdata2)
ydata.mean <- getdum$ydata.mean
ydata.variance <- getdum$ydata.variance
cdata2 <- getdum$cdata
Yvar <- colnames(ydata.mean)[-1]
Cvar <- colnames(cdata2)[-1]
bvar <- all.vars(object$random)
ny <- nrow(ynewdata)
nc <- nrow(cnewdata)
Ny <- nrow(ydata2)
Nc <- nrow(cdata2)
Sig <- object$Sig
p1a <- ncol(Sig) - 1
quadmethod <- object$method
getGH <- GetGHmatrix(quadpoint = quadpoint, p1a = p1a)
xsmatrix <- getGH$xsmatrix
wsmatrix <- getGH$wsmatrix
nsig <- p1a + 1
ynewdata.mean <- ydata.mean[c((Ny-ny+1):Ny), ]
ynewdata.variance <- ydata.variance[c((Ny-ny+1):Ny), ]
cnewdata <- cdata2[c((Nc-nc+1):Nc), ]
if (length(bvar) > 1) bvar1 <- bvar[1:(length(bvar) - 1)]
yID <- unique(ynewdata.mean[, ID])
N.ID <- length(yID)
cID <- cnewdata[, ID]
if (prod(yID == cID) == 0) {
stop("The order of subjects in ydata doesn't match with cnewdata.")
}
if (!is.null(Last.time)) {
if (is.character(Last.time)) {
if (Last.time %in% colnames(cnewdata)) {
Last.time <- cnewdata[, Last.time]
} else {
stop(paste(Last.time, "is not found in cnewdata."))
}
}
if (is.numeric(Last.time) && (length(Last.time) != nrow(cnewdata)))
stop("The last.time vector does not match cnewdata.")
} else {
Last.time <- cnewdata[, Cvar[1]]
}
Pred <- list()
CompetingRisk <- object$CompetingRisk
if (object$CompetingRisk) {
beta <- object$beta
tau <- object$tau
gamma1 <- object$gamma1
gamma2 <- object$gamma2
alpha1 <- object$alpha1
alpha2 <- object$alpha2
nu1 <- object$vee1
nu2 <- object$vee2
H01 <- object$H01
H02 <- object$H02
Sig <- object$Sig
Predraw1 <- matrix(0, nrow = nrow(cnewdata), ncol = length(u))
Predraw2 <- matrix(0, nrow = nrow(cnewdata), ncol = length(u))
y.obs <- list()
lengthu <- length(u)
for (j in 1:N.ID) {
subNDy.mean <- ynewdata.mean[ynewdata.mean[, ID] == yID[j], ]
subNDy.variance <- ynewdata.variance[ynewdata.variance[, ID] == yID[j], ]
subNDc <- cnewdata[cnewdata[, ID] == yID[j], ]
y.obs[[j]] <- data.frame(ynewdata[ynewdata[, ID] == yID[j], c(obs.time, Yvar[1])])
s <- as.numeric(Last.time[j])
CH01 <- CH(H01, s)
CH02 <- CH(H02, s)
Y <- subNDy.mean[, Yvar[1]]
X <- subNDy.mean[, Yvar[2:length(Yvar)]]
X <- as.matrix(X)
W <- subNDy.variance[, -1]
W <- as.matrix(W)
if (nsig == 2) {
Z <- matrix(1, ncol = 1, nrow = length(Y))
} else {
Z <- data.frame(1, subNDy.mean[, bvar1])
Z <- as.matrix(Z)
}
X2 <- as.matrix(subNDc[1, Cvar[3:length(Cvar)]])
if (method == "GH") {
for (jj in 1:lengthu) {
## calculate the CIF
if (quadmethod == "standard") {
CIF <- getECIF(beta, tau, gamma1, gamma2, alpha1, alpha2, nu1,
nu2, Sig, Z, X, W, Y, as.vector(X2), H01, H02,
xsmatrix, wsmatrix, CH01, CH02, s, u[jj])
} else {
data <- list(Y, X, Z, W, X2, CH01, CH02, beta, tau, gamma1, gamma2, alpha1, alpha2, nu1, nu2, Sig)
names(data) <- c("Y", "X", "Z", "W", "X2", "CH01", "CH02", "beta", "tau",
"gamma1", "gamma2", "alpha1", "alpha2", "nu1", "nu2", "Sig")
opt <- optim(rep(0, nsig), logLikCR, data = data, method = "BFGS", hessian = TRUE)
Posmean <- opt$par
PosCov <- solve(opt$hessian)
CIF <- getECIFad(beta, tau, gamma1, gamma2, alpha1, alpha2, nu1,
nu2, Sig, Z, X, W, Y, as.vector(X2), H01, H02,
xsmatrix, wsmatrix, CH01, CH02, s, u[jj], Posmean, PosCov)
}
P1us <- CIF$CIF1
P2us <- CIF$CIF2
Predraw1[j, jj] <- P1us
Predraw2[j, jj] <- P2us
}
} else {
data <- list(Y, X, Z, W, X2, CH01, CH02, beta, tau, gamma1, gamma2, alpha1, alpha2, nu1, nu2, Sig)
names(data) <- c("Y", "X", "Z", "W", "X2", "CH01", "CH02", "beta", "tau",
"gamma1", "gamma2", "alpha1", "alpha2", "nu1", "nu2", "Sig")
opt <- optim(rep(0, nsig), logLikCR, data = data, method = "BFGS", hessian = TRUE)
meanbw <- opt$par
for (jj in 1:lengthu) {
## calculate the CIF
CIF1 <- CIF1.CR(data, H01, H02, s, u[jj], meanbw)
P1us <- Pk.us(CIF1, data, meanbw)
Predraw1[j, jj] <- P1us
CIF2 <- CIF2.CR(data, H01, H02, s, u[jj], meanbw)
P2us <- Pk.us(CIF2, data, meanbw)
Predraw2[j, jj] <- P2us
}
quadpoint = NULL
}
}
for (jj in 1:N.ID) {
Pred[[jj]] <- data.frame(u, Predraw1[jj, ], Predraw2[jj, ])
colnames(Pred[[jj]]) <- c("times", "CIF1", "CIF2")
}
} else {
Predraw <- matrix(0, nrow = nrow(cnewdata), ncol = length(u))
beta <- object$beta
tau <- object$tau
gamma <- object$gamma1
alpha <- object$alpha1
nu <- object$vee1
H01 <- object$H01
Sig <- object$Sig
y.obs <- list()
lengthu <- length(u)
for (j in 1:N.ID) {
subNDy.mean <- ynewdata.mean[ynewdata.mean[, ID] == yID[j], ]
subNDy.variance <- ynewdata.variance[ynewdata.variance[, ID] == yID[j], ]
subNDc <- cnewdata[cnewdata[, ID] == yID[j], ]
y.obs[[j]] <- data.frame(ynewdata[ynewdata[, ID] == yID[j], c(obs.time, Yvar[1])])
CH0 <- CH(H01, Last.time[j])
CH0u <- vector()
for (jj in 1:lengthu) {
CH0u[jj] <- CH(H01, u[jj])
}
Y <- subNDy.mean[, Yvar[1]]
X <- subNDy.mean[, Yvar[2:length(Yvar)]]
X <- as.matrix(X)
W <- subNDy.variance[, -1]
W <- as.matrix(W)
if (nsig == 2) {
Z <- matrix(1, ncol = 1, nrow = length(Y))
} else {
Z <- data.frame(1, subNDy.mean[, bvar1])
Z <- as.matrix(Z)
}
X2 <- as.matrix(subNDc[1, Cvar[3:length(Cvar)]])
if (method == "Laplace") {
## find out E(theta_i)
data <- list(Y, X, Z, W, X2, CH0, beta, tau, gamma, alpha, nu, Sig)
names(data) <- c("Y", "X", "Z", "W", "X2", "CH0", "beta", "tau", "gamma", "alpha", "nu", "Sig")
opt <- optim(rep(0, nsig), logLik, data = data, method = "BFGS", hessian = TRUE)
meanbw <- opt$par
for (jj in 1:lengthu) {
Pi <- P.us(data, CH0u[jj], meanbw)
Predraw[j, jj] <- 1 - Pi
}
quadpoint <- NULL
} else {
for (jj in 1:lengthu) {
if (quadmethod == "standard") {
Predraw[j, jj] <- getES(beta, tau, gamma, alpha, nu, Sig, Z, X, W, Y,
as.vector(X2), xsmatrix, wsmatrix, CH0, CH0u[jj])
} else {
data <- list(Y, X, Z, W, X2, CH0, beta, tau, gamma, alpha, nu, Sig)
names(data) <- c("Y", "X", "Z", "W", "X2", "CH0", "beta", "tau",
"gamma1", "alpha1", "nu1", "Sig")
opt <- optim(rep(0, nsig), logLik, data = data, method = "BFGS", hessian = TRUE)
Posmean <- opt$par
PosCov <- solve(opt$hessian)
Predraw[j, jj] <- getESad(beta, tau, gamma, alpha, nu, Sig, Z, X, W, Y,
as.vector(X2), xsmatrix, wsmatrix, CH0, CH0u[jj], Posmean, PosCov)
}
}
}
}
for (jj in 1:N.ID) {
Pred[[jj]] <- data.frame(u, Predraw[jj, ])
colnames(Pred[[jj]]) <- c("times", "PredSurv")
}
}
names(y.obs) <- names(Pred) <- yID
Last.time <- data.frame(cID, Last.time)
colnames(Last.time)[1] <- ID
sum <- list(Pred = Pred, Last.time = Last.time, y.obs = y.obs, method = method, quadpoint = quadpoint,
CompetingRisk = CompetingRisk, quadmethod = quadmethod)
class(sum) <- "survfitJMMLSM"
sum
}
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.