Nothing
#' Simulate data from the Fine-Gray Model
#'
#' @description Simulate data from the model proposed in Fine and Gray (1999) for two causes. Cause 1 is assumed
#' to be of primary importance.
#'
#' @param nobs Integer: Number of observations in simulated dataset.
#' @param beta1 A vector of effect sizes for cause 1 of length ncovs
#' @param beta2 A vector of effect sizes for cause 2 of length ncovs
#' @param X A matrix of fixed covariates (nobs x ncovs). If \code{X} is NULL (default) then \code{X} will be simulated from MVN(O, I) with n = nobs and p = length(beta1).
#' @param u.min Numeric: controls lower bound of censoring distribution where C ~ U(u.min, u.max)
#' @param u.max Numeric: controls upper bound of censoring distribution where C ~ U(u.min, u.max)
#' @param p Numeric: value between 0 and 1 which controls the mixture probability.
#' @param returnX Logical: Whether to return \code{X} or not. Default is TRUE. Recommended if \code{X} is NULL.
#' @details The function simulates data according to the setup by Fine and Gray (1999). See their paper for more information.
#' @return Returns a list with the following:
#' \item{ftime}{vector of \code{nobs} simulated event times}
#' \item{ftime}{vector of \code{nobs} simulated event indicators (0/1/2)}
#' \item{X}{design matrix if \code{returnX = TRUE}. (simulated design matrix if \code{X = NULL}.)}
#' @import stats
#' @examples
#' set.seed(2019)
#' nobs <- 500
#' beta1 <- c(0.40, -0.40, 0, -0.50, 0, 0.60, 0.75, 0, 0, -0.80)
#' beta2 <- -beta1
#' Z <- matrix(rnorm(nobs * length(beta1)), nrow = nobs)
#' dat <- simulateTwoCauseFineGrayModel(nobs, beta1, beta2, Z, u.min = 0, u.max = 1, p = 0.5)
#' @references
#' Fine J. and Gray R. (1999) A proportional hazards model for the subdistribution of a competing risk. \emph{JASA} 94:496-509.
#' @export
simulateTwoCauseFineGrayModel <- function(nobs, beta1, beta2, X = NULL, u.min = 0, u.max, p = 0.5,
returnX = FALSE) {
#warnings:
if(length(beta1) != length(beta2)) stop("Dimension of beta1 and beta2 should be the same")
ncovs <- length(beta1)
if(is.null(X)) {
X <- matrix(rnorm(nobs * ncovs), nrow = nobs)
returnX <- TRUE
}
#- Generate indicator for cause
c.ind <- 1 + rbinom(nobs, 1, prob = (1 - p)^exp(X %*% beta1))
#Conditional on cause indicators, we simulate the model.
ftime <- numeric(nobs)
eta1 <- X[c.ind == 1, ] %*% beta1 #linear predictor for cause on interest
eta2 <- X[c.ind == 2, ] %*% beta2 #linear predictor for competing risk
u1 <- runif(length(eta1))
t1 <- -log(1 - (1 - (1 - u1 * (1 - (1 - p)^exp(eta1)))^(1 / exp(eta1))) / p)
t2 <- rexp(length(eta2), rate = exp(eta2))
ci <- runif(nobs, min = u.min, max = u.max) #simulate censoring times
ftime[c.ind == 1] <- t1
ftime[c.ind == 2] <- t2
ftime <- pmin(ftime, ci)
fstatus <- ifelse(ftime == ci, 0, 1)
fstatus <- fstatus * c.ind
out <- list()
out$ftime <- ftime
out$fstatus <- fstatus
if(returnX) out$X <- X
return(out)
}
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.