Nothing
R/sim-pc-hazard.R
In mets: Analysis of Multivariate Event Times
Defines functions
read.fit
pre.cifs
sim.cifsRestrict
sim.cifs
sim.cif
subdist
invsubdist
simsubdist
sim.cause.cox
sim.base
sim.cox
cause.pchazard.sim
rcrisks
rcrisk
simrchaz
addCums
lin.approx
rchaz
Documented in
addCums
cause.pchazard.sim
invsubdist
lin.approx
pre.cifs
rchaz
rcrisk
rcrisks
read.fit
sim.base
sim.cause.cox
sim.cif
sim.cifs
sim.cifsRestrict
sim.cox
simrchaz
simsubdist
subdist
#' Simulation of Piecewise constant hazard model (Cox).
#'
#' Simulates data from piecwise constant baseline hazard that can also be of
#' Cox type. Censor data at highest value of the break points.
#'
#' For a piecewise linear cumulative hazard the inverse is easy to compute with
#' and delayed entry x we compute
#' \deqn{\Lambda^{-1}(\Lambda(x) + E/RR)},
#' where RR are the relative risks and E is exponential with mean 1.
#' This quantity has survival function
#' \deqn{P(T > t | T>x) = exp(-RR (\Lambda(t) - \Lambda(x)))}.
#'
#' @param cumhazard cumulative hazard, or piece-constant rates for periods defined by first column of input.
#' @param rr relative risk for simulations, alternatively when rr=1 specify n
#' @param n number of simulation if rr not given
#' @param entry delayed entry time for simuations.
#' @param cum.hazard specifies wheter input is cumulative hazard or rates.
#' @param cause name of cause
#' @param extend to extend piecewise constant with constant rate. Default is average rate over time from cumulative (when TRUE), if numeric then uses given rate.
#' @author Thomas Scheike
#' @keywords survival
#' @examples
#'
#' chaz <- c(0,1,1.5,2,2.1)
#' breaks <- c(0,10, 20, 30, 40)
#' cumhaz <- cbind(breaks,chaz)
#' n <- 100
#' X <- rbinom(n,1,0.5)
#' beta <- 0.2
#' rrcox <- exp(X * beta)
#'
#' pctime <- rchaz(cumhaz,n=1000,cum.hazard=FALSE)
#' pctimecox <- rchaz(cumhaz,rrcox,cum.hazard=FALSE)
#'
#' @export
#' @aliases simrchaz addCums lin.approx
rchaz <- function(cumhazard,rr,n=NULL,entry=NULL,cum.hazard=TRUE,cause=1,extend=FALSE)
{# {{{
### cumh=cbind(breaks,rates), first rate is 0 if cumh=FALSE
### cumh=cbind(breaks,cumhazard) if cumh=TRUE
### if (is.null(rr) & is.null(n)) stop("must give rr or n\n");
if (!is.null(n)) rr <- rep(1,n)
n <- length(rr)
breaks <- cumhazard[,1]
rates <- cumhazard[,2][-1]
mm <- tail(breaks,1)
if (cum.hazard==FALSE) {
cumh <- cumsum(c(0,diff(breaks)*rates))
cumhazard <- cbind(breaks,cumh)
} else cumh <- cumhazard[,2]
ttt <- rexp(n)/rr
if (cumhazard[1,2]>0) { ## start cumulative hazard with a 0
### warning("Safest to start with cumulative hazard 0 to avoid problems\n");
cumhazard <- rbind(c(0,0),cumhazard)
cumh <- c(0,cumh)
}
###
if (!is.null(entry)) {
if (length(entry)==1) entry <- rep(entry,n) else entry <- entry
cumentry <- lin.approx(entry,cumhazard,x=1)
} else { entry <- cumentry <- rep(0,n) }
###
ttte <- ttt+cumentry
rrx <- lin.approx(ttte,cumhazard,x=-1)
rrx <- ifelse(rrx>mm,mm,rrx)
status <- rep(0,n)
status <- ifelse(rrx<mm,cause,status)
tcum <- tail(cumhazard,1)
extend.rate <- NULL
if (is.logical(extend))
if (extend) extend.rate <- tcum[2]/tcum[1]
if (is.numeric(extend)) { extend.rate <- extend; extend <- TRUE;}
if (extend) {
whoe <- which(status==0)
nn <- length(whoe)
if (nn>0) {
textend <- rexp(nn)/(rr[whoe]*extend.rate)
rrx[whoe] <- mm+textend
status[whoe] <- 1
}
}
dt <- data.frame(entry=entry,time=rrx,status=status,rr=rr)
colnames(dt) <- c("entry","time","status","rr");
attr(dt,"cumhaz") <- cumhazard
attr(dt,"extend.rate") <- extend.rate
return(dt)
}# }}}
#' @export
lin.approx <- function(x2,xfx,x=1)
{# {{{
### x=1 gives f(x2)
### x=-1 gives f^-1(x2)
breaks <- xfx[,x]
fx <- xfx[,-x]
ri <- fast.approx(breaks,x2,type="left")
maxindex <- which(ri==length(breaks))
rip1 <- ri+1
rip1[maxindex] <- length(breaks)
rrr <- (x2-breaks[ri])/(breaks[rip1]-breaks[ri])
rrr[maxindex] <- 0
res <- rrr*(fx[rip1]-fx[ri])+fx[ri]
res[is.na(res)] <- tail(fx,1)
return(res)
}# }}}
#' @export
addCums <- function(cumB,cumA,max=NULL)
{# {{{
## take times
times <- sort(unique(c(cumB[,1],cumA[,1])))
if (!is.null(max)) times <- times[times<max]
cumBjx <- lin.approx(times,cumB,x=1)
cumAjx <- lin.approx(times,cumA,x=1)
cumBA <- cumBjx+cumAjx
return(cbind(times,cumBA))
}# }}}
#' @export
simrchaz <- function(cumhazard,rr,n=NULL,cens=NULL,rrc=NULL,...)
{# {{{
### adds censoring to to rchaz call
if (!is.null(n)) rr <- rep(1,n)
n <- length(rr)
ptt <- rchaz(cumhazard,rr,...)
if (is.null(rrc)) {
rrc <- rep(1,n)
}
if (is.matrix(cens)) {
pct <- rchaz(cens,rrc,...)
pct <- pct$time
} else {
if (is.numeric(cens)) pct<- rexp(n)/cens else {
chaz <-sum(ptt$status)/sum(ptt$time) ## hazard averate T haz
pct<- rexp(n)/chaz
}
}
dt <- data.frame(time=pmin(ptt$time,pct), status=ifelse(ptt$time<pct,ptt$status,0))
return(dt)
}# }}}
#' Simulation of Piecewise constant hazard models with two causes (Cox).
#'
#' Simulates data from piecwise constant baseline hazard that can also be of
#' Cox type. Censor data at highest value of the break points for either of the
#' cumulatives.
#'
#' @param cumhaz1 cumulative hazard of cause 1
#' @param cumhaz2 cumulative hazard of cause 1
#' @param rr1 number of simulations or vector of relative risk for simuations.
#' @param rr2 number of simulations or vector of relative risk for simuations.
#' @param n number of simulation if rr not given
#' @param cens to censor further , rate or cumumlative hazard
#' @param rrc retlativ risk for censoring.
#' @param ... arguments for rchaz
#' @author Thomas Scheike
#' @keywords survival
#' @examples
#' library(mets); data(bmt); library(survival)
#'
#' cox1 <- phreg(Surv(time,cause==1)~tcell+platelet,data=bmt)
#' cox2 <- phreg(Surv(time,cause==2)~tcell+platelet,data=bmt)
#'
#' X1 <- bmt[,c("tcell","platelet")]
#' n <- 100
#' xid <- sample(1:nrow(X1),n,replace=TRUE)
#' Z1 <- X1[xid,]
#' Z2 <- X1[xid,]
#' rr1 <- exp(as.matrix(Z1) %*% cox1$coef)
#' rr2 <- exp(as.matrix(Z2) %*% cox2$coef)
#'
#' d <- rcrisk(cox1$cum,cox2$cum,rr1,rr2)
#' dd <- cbind(d,Z1)
#'
#' scox1 <- phreg(Surv(time,status==1)~tcell+platelet,data=dd)
#' scox2 <- phreg(Surv(time,status==2)~tcell+platelet,data=dd)
#' par(mfrow=c(1,2))
#' plot(cox1); plot(scox1,add=TRUE)
#' plot(cox2); plot(scox2,add=TRUE)
#' cbind(cox1$coef,scox1$coef,cox2$coef,scox2$coef)
#'
#' @export
#' @aliases cause.pchazard.sim rcrisks
rcrisk <-function(cumhaz1,cumhaz2,rr1,rr2,n=NULL,cens=NULL,rrc=NULL,...)
{#'# {{{
if (!is.null(n)) { rr1 <- rep(1,n); rr2 <- rep(1,n) }
n <- length(rr1);
if (missing(rr2)) rr2 <-rep(1,n)
ptt1 <- rchaz(cumhaz1,rr1,cum.hazard=TRUE,...)
ptt2 <- rchaz(cumhaz2,rr2,cum.hazard=TRUE,...)
ptt <- data.frame(time=pmin(ptt1$time,ptt2$time),status=ifelse(ptt1$time<=ptt2$time,ptt1$status,ptt2$status*2),entry=ptt1$entry)
if (!is.null(cens)) {
if (is.null(rrc)) rrc <- rep(1,n)
if (is.matrix(cens)) {
pct <- rchaz(cens,rrc,...)
pct <- pct$time
} else {
if (is.numeric(cens)) pct<- rexp(n)/cens else {
chaz <-sum(ptt1$status+ptt2$status)/sum(ptt1$time+ptt2$time) ## hazard averate T haz
pct<- rexp(n)/chaz
}
}
ptt <- data.frame(time=pmin(ptt$time,pct),status=ifelse(ptt$time<pct,ptt$status,0),entry=ptt$entry)
}
return(ptt)
}# }}}
#' @export
rcrisks <-function(cumhazs,rrs,n=NULL,cens=NULL,rrc=NULL,entry=NULL,causes=NULL,...)
{#'# {{{
status <- NULL
if (!is.null(n)) rrs <- matrix(1,n,length(cumhazs))
n <- nrow(rrs);
if (!is.null(cens)) cens <- list(cens)
if (!is.null(cens)) {
if (is.null(rrc)) rrc <- rep(1,n);
rrs <- cbind(rrs,rrc)
}
cumss <- c(cumhazs,cens)
for (i in seq_along(cumss)) cumss[[i]] <- rbind(0,cumss[[i]])
cumss <- extendCums(cumss,NULL)
## first time
ptt <- rchaz(cumss[[1]],rrs[,1],entry=entry)
## other times and always min
if (length(cumss)>=2)
for (i in 2:length(cumss)) {
cum <- cumss[[i]]
ptt1 <- rchaz(cumss[[i]],rrs[,i],entry=entry)
ptt <- data.frame(time=pmin(ptt$time,ptt1$time),entry=ptt1$entry,
status=ifelse(ptt$time<=ptt1$time,ptt$status,ptt1$status*i))
}
if (!is.null(cens)) ptt <- dtransform(ptt,status=0,status==length(cumss))
if (!is.null(causes)) {
ptt$status <- c(0,causes)[ptt$status+1]
}
return(ptt)
}# }}}
#' @export
cause.pchazard.sim<-function(cumhaz1,cumhaz2,rr1,rr2,cens=NULL,rrc=NULL,...)
{
rcrisk(cumhaz1,cumhaz2,rr1,rr2,cens=NULL,rrc=NULL,...)
}
#' Simulation of output from Cox model.
#'
#' Simulates data that looks like fit from Cox model. Censor data automatically
#' for highest value of the event times by using cumulative hazard.
#'
#' @param cox output form coxph or cox.aalen model fitting cox model.
#' @param n number of simulations.
#' @param data to extract covariates for simulations (draws from observed
#' covariates).
#' @param cens specifies censoring model, if "is.matrix" then uses cumulative
#' hazard given, if "is.scalar" then uses rate for exponential, and if not
#' given then takes average rate of in simulated data from cox model.
#' @param rrc possible vector of relative risk for cox-type censoring.
#' @param entry delayed entry variable for simulation.
#' @param ... arguments for rchaz, for example entry-time
#' @author Thomas Scheike
#' @keywords survival
#' @examples
#'
#' data(sTRACE)
#' cox <- coxph(Surv(time,status==9)~vf+chf+wmi,data=sTRACE)
#' sim1 <- sim.cox(cox,1000,data=sTRACE)
#' cc <- coxph(Surv(time,status)~vf+chf+wmi,data=sim1)
#' cbind(cox$coef,cc$coef)
#'
#' cor(sim1[,c("vf","chf","wmi")])
#' cor(sTRACE[,c("vf","chf","wmi")])
#'
#' cox <- phreg(Surv(time, status==9)~vf+chf+wmi,data=sTRACE)
#' sim3 <- sim.cox(cox,1000,data=sTRACE)
#' cc <- phreg(Surv(time, status)~vf+chf+wmi,data=sim3)
#' cbind(cox$coef,cc$coef)
#' plot(cox,se=TRUE)
#' plot(cc,add=TRUE,col=2)
#'
#' cox <- phreg(Surv(time,status==9)~strata(chf)+vf+wmi,data=sTRACE)
#' sim3 <- sim.cox(cox,100,data=sTRACE)
#' cc <- phreg(Surv(time, status)~strata(chf)+vf+wmi,data=sim3)
#' cbind(cox$coef,cc$coef)
#' plot(cox)
#' plot(cc,add=TRUE,col=2)
#'
#' @aliases read.fit sim.base simulate.cox
#' @export sim.cox
#' @usage sim.cox(cox,n,data=NULL,cens=NULL,rrc=NULL,entry=NULL,...)
sim.cox <- function(cox,n,data=NULL,cens=NULL,rrc=NULL,entry=NULL,...)
{# {{{
### cumh=cbind(breaks,rates), first rate is 0 if cumh=FALSE
### cumh=cbind(breaks,cumhazard) if cumh=TRUE
des <- read.fit(cox,n,data=data,...)
Z <- des$Z; cumhazard <- des$cum; rr <- des$rr;
if (!inherits(cox,"phreg")) {
if (cumhazard[1,2]>0) cumhazard <- rbind(c(0,0),cumhazard)
ptt <- rchaz(cumhazard,rr,entry=entry)
ptt <- cbind(ptt,Z)
} else {
ptt <- data.frame()
stratj <- cox$strata[cox$jumps]
if (cox$nstrata>1) {
for (j in 0:(cox$nstrata-1)) {
cumhazardj <- rbind(c(0,0),cox$cumhaz[stratj==j,])
if (!is.null(entry)) entry <- entry[des$strataid==j]
pttj <- rchaz(cumhazardj,rr[des$strataid==j],entry=entry)
Zj <- Z[des$strataid==j,,drop=FALSE]
pttj <- cbind(pttj,Zj)
ptt <- rbind(ptt,pttj)
}
} else {
if (cumhazard[1,2]>0) cumhazard <- rbind(c(0,0),cumhazard)
ptt <- rchaz(cumhazard,rr,entry=entry)
ptt <- cbind(ptt,Z)
}
}
if (!is.null(cens)) {
if (is.null(rrc)) rrc <- rep(1,n)
if (is.matrix(cens)) {
pct <- rchaz(cens,rrc,entry=entry)
pct <- pct$time
}
else {
if (is.numeric(cens)) pct<- rexp(n)/cens else {
chaz <-sum(ptt$status)/sum(ptt$time) ## hazard averate T haz
pct<- rexp(n)/chaz
if (!is.null(entry)) pct <- entry + pct
}
}
ptt$time <- pmin(ptt$time,pct)
ptt$status <- ifelse(ptt$time<pct,ptt$status,0)
}
attr(ptt,"id") <- des$id
return(ptt)
}# }}}
#' @export sim.base
#' @usage sim.base(cumhaz,n,data=NULL,strata=NULL,cens=NULL,entry=NULL,...)
sim.base <- function(cumhaz,n,data=NULL,strata=NULL,cens=NULL,entry=NULL,...)
{# {{{
### cumh=cbind(breaks,rates), first rate is 0 if cumh=FALSE
### cumh=cbind(breaks,cumhazard) if cumh=TRUE
if (!is.null(strata)) {
us <- unique(strata)
nus <- length(us)
if (length(n)!=nus) warning("n must be given for each strata as vector \n")
if (length(n)!=nus) n <- rep(round(n/nus),nus)
}
nt <- 0
if (is.null(strata)) {
if (cumhaz[1,2]>0) cumhaz <- rbind(c(0,0),cumhaz)
ptt <- rchaz(cumhaz,n=n,entry=entry)
nt <- n
} else {
ptt <- c()
nstrata <- length(unique(strata))
for (i in seq_along(unique(strata))) {
j <- unique(strata)[i]
cumhazardj <- rbind(c(0,0),cumhaz[strata==j,])
if (!is.null(entry)) entry <- entry[strata==j]
ns <- n[i]
nt <- nt+ns
pttj <- cbind(rchaz(cumhazardj,n=ns,entry=entry),j)
colnames(pttj)[5] <- "strata"
ptt <- rbind(ptt,pttj)
}
}
if (!is.null(cens)) {
if (is.matrix(cens)) {
pct <- rchaz(cens,n=nt,entry=entry)$time
}
else {
pct<- rexp(nt)/(cens)
}
ptt$time <- pmin(ptt$time,pct)
ptt$status <- ifelse(ptt$time<pct,ptt$status,0)
}
return(ptt)
}# }}}
#' Simulation of cause specific from Cox models.
#'
#' Simulates data that looks like fit from cause specific Cox models.
#' Censor data automatically. When censoring is given in the list of causes this
#' will give censoring that looks like the data. Covariates are drawn from data-set
#' with replacement. This gives covariates like the data.
#'
#'
#' @param coxs list of cox models.
#' @param n number of simulations.
#' @param data to extract covariates for simulations (draws from observed
#' covariates).
#' @param cens specifies censoring model, if NULL then only censoring for
#' each cause at end of last event of this type.
#' if "is.matrix" then uses cumulative.
#' hazard given, if "is.scalar" then uses rate for exponential, and if not
#' given then takes average rate of in simulated data from cox model.
#' But censoring can also be given as a cause.
#' @param rrc possible vector of relative risk for cox-type censoring.
#' @param ... arguments for rchaz, for example entry-time
#' @author Thomas Scheike
#' @keywords survival
#' @examples
#' nsim <- 100
#'
#' data(bmt)
#' # coxph
#' cox1 <- coxph(Surv(time,cause==1)~tcell+platelet,data=bmt)
#' cox2 <- coxph(Surv(time,cause==2)~tcell+platelet,data=bmt)
#' coxs <- list(cox1,cox2)
#' dd <- sim.cause.cox(coxs,nsim,data=bmt)
#' scox1 <- coxph(Surv(time,status==1)~tcell+platelet,data=dd)
#' scox2 <- coxph(Surv(time,status==2)~tcell+platelet,data=dd)
#' cbind(cox1$coef,scox1$coef)
#' cbind(cox2$coef,scox2$coef)
#'
#' data(bmt)
#' cox1 <- phreg(Surv(time,cause==1)~tcell+platelet,data=bmt)
#' cox2 <- phreg(Surv(time,cause==2)~tcell+platelet,data=bmt)
#' coxs <- list(cox1,cox2)
#' dd <- sim.cause.cox(coxs,nsim,data=bmt)
#' scox1 <- phreg(Surv(time,status==1)~tcell+platelet,data=dd)
#' scox2 <- phreg(Surv(time,status==2)~tcell+platelet,data=dd)
#' cbind(cox1$coef,scox1$coef)
#' cbind(cox2$coef,scox2$coef)
#' par(mfrow=c(1,2))
#' plot(cox1); plot(scox1,add=TRUE);
#' plot(cox2); plot(scox2,add=TRUE);
#'
#' cox1 <- phreg(Surv(time,cause==1)~strata(tcell)+platelet,data=bmt)
#' cox2 <- phreg(Surv(time,cause==2)~strata(tcell)+platelet,data=bmt)
#' coxs <- list(cox1,cox2)
#' dd <- sim.cause.cox(coxs,nsim,data=bmt)
#' scox1 <- phreg(Surv(time,status==1)~strata(tcell)+platelet,data=dd)
#' scox2 <- phreg(Surv(time,status==2)~strata(tcell)+platelet,data=dd)
#' cbind(cox1$coef,scox1$coef)
#' cbind(cox2$coef,scox2$coef)
#' par(mfrow=c(1,2))
#' plot(cox1); plot(scox1,add=TRUE);
#' plot(cox2); plot(scox2,add=TRUE);
#'
#' @export sim.cause.cox
#' @usage sim.cause.cox(coxs,n,data=NULL,cens=NULL,rrc=NULL,...)
sim.cause.cox <- function(coxs,n,data=NULL,cens=NULL,rrc=NULL,...)
{# {{{
### cumh=cbind(breaks,rates), first rate is 0 if cumh=FALSE
### cumh=cbind(breaks,cumhazard) if cumh=TRUE
if (!is.list(coxs)) stop("Cox models in list form\n");
ptt <- sim.cox(coxs[[1]],n,data=data)
simcovs <- ptt[,(5:ncol(ptt))]
i=2
if (length(coxs)>=2)
for (i in 2:length(coxs)) {
cox <- coxs[[i]]
ptt1 <- sim.cox(coxs[[i]],n,data=data,id=attr(ptt,"id"))
Zi <- ptt1[,(5:ncol(ptt1))]
pm <- match(names(Zi),names(simcovs))
Zi <- Zi[,-pm]
simcovs <- cbind(simcovs,Zi)
dt <- data.frame(time=pmin(ptt$time,ptt1$time),
status=ifelse(ptt$time<=ptt1$time,ptt$status,ptt1$status*i))
}
dt <- cbind(dt,simcovs)
if (!is.null(cens)) {
if (is.matrix(cens)) {
pct <- rchaz(cens,rrc)
pct <- pct$time
}
else {
if (is.numeric(cens)) pct<- rexp(n)/cens else {
chaz <-sum(ptt$status)/sum(ptt$time) ## hazard averate T haz
pct<- rexp(n)/chaz
}
}
dt <- cbind(data.frame(time=pmin(ptt$time,pct),status=ifelse(ptt$time<pct,ptt$status,0)),simcovs)
}
return(dt)
}# }}}
#' @export
simsubdist <- function(cumhazard,rr,n=NULL,entry=NULL,type="cloglog",startcum=c(0,0),...)
{# {{{
## Fine-Gray model cloglog F1= 1-exp(-cum(t)*rr)
## logistic F1= cum(t)*rr/(1+cum(t)*rr)
## identity F1= cum(t)*rr
## rr=exp(X^t beta)
if (!is.null(n)) rr <- rep(1,n)
entry=NULL
logit <- function(p) log(p/(1-p))
breaks <- cumhazard[,1]
rates <- cumhazard[,2][-1]
mm <- tail(breaks,1)
cumh <- cumhazard[,2]
n <- length(rr)
if (cumh[1]>0) {
cumh <- c(startcum[2],cumh)
breaks <- c(startcum[1],breaks)
}
if (type=="cloglog") {
F1tau <- 1-exp(-tail(cumh,1)*rr)
ttt <- -log(1-runif(n)*F1tau)/rr
} else if (type=="logistic") {
F1tau <- tail(cumh,1)*rr/(1+tail(cumh,1)*rr)
v <- runif(n)*F1tau
ttt <- exp(logit(v))/rr;
} else if (type=="identity" | type=="cif") {
F1tau <- tail(cumh,1)
ttt <- runif(n)*F1tau
## rr only affects binomial draw
F1tau <- F1tau*rr
} else stop(" cloglog or logistic or give function (fun=) \n");
###
entry <- cumentry <- rep(0,n)
ttte <- ttt+cumentry
rrx <- lin.approx(ttt,cbind(breaks,cumh),x=-1)
timecause <- rrx
###
rrx <- ifelse(rrx>mm,mm,rrx)
status <- rbinom(n,1,F1tau)
rrx[status==0] <- mm
dt <- data.frame(entry=entry,time=rrx,status=status,rr=rr,F1tau=F1tau,timecause=timecause)
attr(dt,"cumhaz") <- cumhazard
return(dt)
}# }}}
#' @export
invsubdist <- function(F1,u,entry=NULL,cond=1,ptrunc=NULL)
{# {{{
### computes inverse subdist F1^{-1}(u)
### computes inverse subdist F1^{-1}(u)
F1max <- tail(F1[,2],1)
mmax <- tail(F1[,1],1)
if (!is.null(entry)) F1entry <- subdist(F1,entry)[,2];
if (!is.null(entry)) if (is.null(ptrunc)) ptrunc <- 1- F1entry
if (cond==0 & !is.null(entry)) u<-u*ptrunc+F1entry
if (cond==1 & !is.null(entry)) u<-u*(F1max-F1entry)+F1entry
if (cond==1 & is.null(entry)) u<-u*F1max
rrx <- lin.approx(u,F1,x=-1)
status <- rep(1,length(u))
status[u>=F1max] <- 0
dt <- data.frame(time=rrx,status=status,y=u)
if (!is.null(entry)) dt <- cbind(dt,entry)
attr(dt,"F1") <- F1
attr(dt,"Fmax") <- F1max
return(dt)
}# }}}
#' @export
subdist <- function(F1,times)
{# {{{
rrx <- lin.approx(times,F1,x=1)
dt <- cbind(times,rrx)
colnames(dt) <- c("time","subdist")
return(dt)
}# }}}
#' Simulation of output from Cumulative incidence regression model
#'
#' Simulates data that looks like fit from fitted cumulative incidence model
#'
#' @param cif output form prop.odds.subdist or ccr (cmprsk), can also call invsubdist with
#' with cumulative and linear predictor
#' @param n number of simulations.
#' @param data to extract covariates for simulations (draws from observed
#' covariates).
#' @param Z to use these covariates for simulation rather than drawing new ones.
#' @param drawZ to random sample from Z or not
#' @param cens specifies censoring model, if "is.matrix" then uses cumulative
#' hazard given, if "is.scalar" then uses rate for exponential, and if not
#' given then takes average rate of in simulated data from cox model.
#' @param rrc possible vector of relative risk for cox-type censoring.
#' @param cumstart to start cumulatives at time 0 in 0.
#' @param ... arguments for invsubdist
#' @author Thomas Scheike
#' @keywords survival
#' @examples
#' data(bmt)
#'
#' scif <- cifreg(Event(time,cause)~tcell+platelet+age,data=bmt,cause=1,prop=NULL)
#' summary(scif)
#' plot(scif)
#' ################################################################
#' # simulating several causes with specific cumulatives
#' ################################################################
#'
#' cif1 <- cifreg(Event(time,cause)~tcell+age,data=bmt,cause=1,prop=NULL)
#' cif2 <- cifreg(Event(time,cause)~tcell+age,data=bmt,cause=2,prop=NULL)
#' # dd <- sim.cifsRestrict(list(cif1,cif2),200,data=bmt)
#
#' dd <- sim.cifs(list(cif1,cif2),200,data=bmt)
#' scif1 <- cifreg(Event(time,cause)~tcell+age,data=dd,cause=1)
#' scif2 <- cifreg(Event(time,cause)~tcell+age,data=dd,cause=2)
#'
#' par(mfrow=c(1,2))
#' plot(cif1); plot(scif1,add=TRUE,col=2)
#' plot(cif2); plot(scif2,add=TRUE,col=2)
#' @export sim.cif
#' @aliases sim.cif sim.cifs subdist pre.cifs sim.cifsRestrict simsubdist invsubdist
#' @usage sim.cif(cif,n,data=NULL,Z=NULL,drawZ=TRUE,cens=NULL,rrc=NULL,cumstart=c(0,0),...)
sim.cif <- function(cif,n,data=NULL,Z=NULL,drawZ=TRUE,cens=NULL,rrc=NULL,cumstart=c(0,0),...)
{# {{{
## also extracts coefficients and baseline from coxph, cox.aalen, phreg
## and uses them assuming that the model is cloglog unless otherwise
if (!inherits(cif,"defined")) {
des <- read.fit(cif,n,data=data,Z=Z,drawZ=drawZ,...)
Z <- des$Z; cumhaz <- des$cum; rr <- des$rr;
id <- des$id
cumhaz <- rbind(cumstart,cumhaz)
znames <- names(Z);
model <- des$model
} else {
cumhaz <- cif$cumhaz
if (is.null(Z)) stop("must give Z")
rr <- exp(as.matrix(Z) %*% cif$coef)
if (is.data.frame(Z)) znames <- names(Z) else znames <- colnames(Z)
model <- cif$model
id <- 1:nrow(Z)
}
if (!inherits(cif,"phreg")) {
if (model=="logistic2" | model=="logistic") ptt <- simsubdist(cumhaz,rr,type="logistic",...) else ptt <- simsubdist(cumhaz,rr,type="cloglog",...)
ptt <- cbind(ptt,Z)
} else { ### phreg class# {{{
ptt <- data.frame()
if (cif$nstrata>1) {
stratj <- cif$strata[cif$jumps]
for (j in 0:(cif$nstrata-1)) {
cumhazardj <- rbind(c(0,0),cif$cumhaz[stratj==j,])
if (model[3])
pttj <- simsubdist(cumhazardj,rr[des$strataid==j],type="cloglog") else
pttj <- simsubdist(cumhazardj,rr[des$strataid==j],type="logistic")
Zj <- Z[des$strataid==j,,drop=FALSE]
pttj <- cbind(pttj,Zj)
ptt <- rbind(ptt,pttj)
ptt <- rbind(ptt,pttj)
}
} else {
if (model[3]) ptt <- simsubdist(cumhaz,rr,type="cloglog") else
ptt <- simsubdist(cumhaz,rr,type="logistic")
ptt <- cbind(ptt,Z)
}
}# }}}
### adds censoring
if (!is.null(cens)) {# {{{
if (is.null(rrc)) rrc <- rep(1,n)
if (is.matrix(cens)) {
pct <- rchaz(cens,rrc,...)
pct <- pct$time
}
else {
if (is.numeric(cens)) pct<- rexp(n)/cens else {
chaz <-sum(ptt$status)/sum(ptt$time) ## hazard averate T haz
pct<- rexp(n)/chaz
}
}
ptt$time <- pmin(ptt$time,pct)
ptt$status <- ifelse(ptt$time<pct,ptt$status,0)
} # }}}
attr(ptt,"model") <- model
attr(ptt,"id") <- id
attr(ptt,"znames") <- znames
return(ptt)
}# }}}
#' @export sim.cifs
sim.cifs <- function(cifs,n,data=NULL,Z=NULL,cens=NULL,rrc=NULL,max.times=NULL,causes=c(1,2),...)
{# {{{
if (!is.list(cifs)) stop("Cif models in list form\n");
## must consider all models out to last observation times or max.times
cifs <- pre.cifs(cifs,max.times=max.times)
tau <- tail(cifs[[1]]$cum[,1],1)
sim1 <- sim.cif(cifs[[1]],n,data=data,Z=Z)
Z <- sim1[,attr(sim1,"znames")]
if (!inherits(cifs[[2]],"defined")) {
sim2p <- read.fit(cifs[[2]],1,data=data,Z=NULL)
Z2 <- data[attr(sim1,"id"),names(sim2p$Z)]
sim2 <- sim.cif(cifs[[2]],n,data=data,Z=Z2,drawZ=FALSE)
} else {
Z2 <- Z[,attr(cifs[[2]],"znames")]
sim2 <- sim.cif(cifs[[2]],n,data=data,Z=Z2)
}
ptot <- sim1$F1tau+sim2$F1tau
###
rt <- rbinom(n,1,pmin(ptot,1))
rb <- rbinom(n,1,sim1$F1tau/ptot)
cause=ifelse(rb==1,1,2)
time=ifelse(cause==causes[1],sim1$timecause,sim2$timecause)
cause <- rt*cause
time[cause==0] <- tau
ptt <- data.frame(time=time,status=cause,cause=cause,ptot=ptot)
Zcovs <- cbind(Z,Z2)
samecovs <- match(names(Z2),names(Z))
Ze <- Zcovs[,-samecovs]
ptt <- cbind(ptt,Ze)
if (!is.null(cens)) {# {{{
if (is.null(rrc)) rrc <- rep(1,n)
if (is.matrix(cens)) {
pct <- rchaz(cens,rrc,...)
pct <- pct$time
}
else {
if (is.numeric(cens)) pct<- rexp(n)/cens else {
chaz <-sum(ptt$status)/sum(ptt$time) ## hazard averate T haz
pct<- rexp(n)/chaz
}
}
ptt$time <- pmin(ptt$time,pct)
ptt$status <- ifelse(ptt$time<pct,ptt$status,0)
} # }}}
return(ptt)
}# }}}
#' @export sim.cifsRestrict
sim.cifsRestrict <- function(cifs,n,data=NULL,Z=NULL,cens=NULL,rrc=NULL,max.times=NULL,causes=c(1,2),...)
{# {{{
if (!is.list(cifs)) stop("Cif models in list form\n");
## must consider all models out to last observation times or max.times
cifs <- pre.cifs(cifs,max.times=max.times)
tau <- tail(cifs[[1]]$cum[,1],1)
sim1 <- sim.cif(cifs[[1]],n,data=data,Z=Z)
Z <- sim1[,attr(sim1,"znames")]
if (!inherits(cifs[[2]],"defined")) {
sim2p <- read.fit(cifs[[2]],1,data=data,Z=NULL)
Z2 <- data[attr(sim1,"id"),names(sim2p$Z)]
sim2 <- sim.cif(cifs[[2]],n,data=data,Z=Z2,drawZ=FALSE)
} else {
Z2 <- Z[,attr(cifs[[2]],"znames")]
sim2 <- sim.cif(cifs[[2]],n,data=data,Z=Z2)
}
ptot <- sim1$F1tau+sim2$F1tau*(1-sim1$F1tau)
###
rt <- rbinom(n,1,pmin(ptot,1))
rb <- rbinom(n,1,sim1$F1tau/ptot)
cause=ifelse(rb==1,1,2)
time=ifelse(cause==causes[1],sim1$timecause,sim2$timecause)
cause <- rt*cause
time[cause==0] <- tau
ptt <- data.frame(time=time,status=cause,cause=cause,ptot=ptot)
ptt <- cbind(ptt,Z)
samecovs <- match(names(Z2),names(Z))
Ze <- Z2[,-samecovs]
ptt <- cbind(ptt,Ze)
if (!is.null(cens)) {# {{{
if (is.null(rrc)) rrc <- rep(1,n)
if (is.matrix(cens)) {
cum.hazard <- TRUE
pct <- rchaz(cens,rrc,...)
pct <- pct$time
}
else {
if (is.numeric(cens)) pct<- rexp(n)/cens else {
chaz <-sum(ptt$status)/sum(ptt$time) ## hazard averate T haz
pct<- rexp(n)/chaz
}
}
ptt$time <- pmin(ptt$time,pct)
ptt$status <- ifelse(ptt$time<pct,ptt$status,0)
} # }}}
return(ptt)
}# }}}
#' @export pre.cifs
pre.cifs <- function(cifs,n,data=NULL,pprint=FALSE,max.times=NULL,...)
{# {{{
if (!is.list(cifs)) stop("Cif models in list form\n");
## iff needed put cumulative baseline in argument cifs$cum
## extracts maximum time of all cumulatives
maxtimes <- rep(0,length(cifs))
for (i in 1:length(cifs)) {
if (inherits(cifs[[i]],"coxph")) {
stop("Use phreg of mets instead\n");
} else if (inherits(cifs[[i]],"crr")) {
cifs[[i]]$cum <- rbind(c(0,0),cbind(cifs[[i]]$uftime,cumsum(cifs[[i]]$bfitj)))
maxtimes[i] <- max(cifs[[i]]$uftime)
} else if (inherits(cifs[[i]],"phreg")) {
cifs[[i]]$cumhaz <- cifs[[i]]$cumhaz
maxtimes[i] <- max(cifs[[i]]$cumhaz[,1])
} else if (inherits(cifs[[i]],"defined")) {
cifs[[i]]$cumhaz <- cifs[[i]]$cumhaz
maxtimes[i] <- max(cifs[[i]]$cumhaz[,1])
} else {
maxtimes[i] <- max(cifs[[i]]$cum[,1])
}
}
mmtimes <- min(maxtimes)
if (is.null(max.times)) mtimes <- mmtimes else mtimes <- max.times
if (mtimes > mmtimes) {
warning("max.time is further out than max for some cumulatives\n");
cat(" Max times for cumulatives in cifs \n");
print(maxtimes)
}
for (i in 1:length(cifs)) {
cums <- cifs[[i]]$cum
keep <- cums[,1]<mtimes
Fmm <- subdist(cums,mtimes)
cums <- cums[keep,,drop=FALSE]
cums <- rbind(cums,Fmm)
cifs[[i]]$cum <- cums
cifs[[i]]$cumhaz <- cums
if (pprint) {
print(head(cums))
print(tail(cums))
}
}
return(cifs)
}# }}}
## reads a coxph, cox.aalen, phreg, crr, comprisk, prop.odds.subdist object
## and returns cumulative hazard, linear predictor of a resample of size
## n of data, for coxph cox.aalen, comprisk the design Z is constructed
## using the data, but Z can also be given which is needed for crr
## if drawZ is true the covariates in Z are used but multplied coefficient
## based on column names, if fixZ=TRUE the matrix Z is multiplied coefficients
## as is Z %*% coef(x) to form linear predictor
read.fit <- function(cox,n,data=NULL,Z=NULL,drawZ=TRUE,fixZ=FALSE,id=NULL)
{# {{{
if (inherits(cox,"coxph"))
{# {{{
## basehaz(cox,centered=FALSE)[,c(2,1)]
sfit <- survival::survfit(cox, se.fit=FALSE)
zcoef <- ifelse(is.na(coef(cox)), 0, coef(cox))
offset <- sum(cox$means * zcoef)
chaz <- sfit$cumhaz * exp(-offset)
base <- cbind(sfit$time,chaz)
mt <- model.frame(cox)
Y <- model.extract(mt, "response")
if (!inherits(Y, "Surv"))
stop("Response must be a survival object")
if (attr(Y, "type") == "right") {
time <- Y[, "time"];
status <- Y[, "status"]
} else stop("Expected right-censored data.");
if (is.null(Z)) Z <- na.omit(model.matrix(cox))
nn <- colnames(Z)
if (fixZ) Z <- Z else Z <- Z[,names(coef(cox)),drop=FALSE]
rownames(Z) <- NULL
jtime <- sort(time[status==1])
cumhazard <- rbind(c(0,0),Cpred(base,jtime))
rr <- exp( Z %*% matrix(coef(cox),ncol=1))
if (drawZ==TRUE) xid <- sample(1:nrow(Z),n,replace=TRUE) else xid <- 1:nrow(Z)
if (!is.null(id)) xid <- id
rr <- rr[xid]
Z <- Z[xid,,drop=FALSE]
colnames(Z) <- nn
model <- "fg"
}# }}}
if (inherits(cox,"cox.aalen"))
{# {{{
formula <- attr(cox, "Formula")
### Z1 <- na.omit(model.matrix(cox,data=data))
p <- nrow(cox$gamma)
if (is.null(Z)) {
Z <- na.omit(get_all_vars(formula,data=data))
nz <- ncol(Z)
Z <- Z[,seq(nz-p+1,nz)]
}
lrr <- as.matrix(Z) %*% cox$gamma
cumhazard <- cox$cum
rr <- exp(lrr)
if (drawZ==TRUE) xid <- sample(1:nrow(Z),n,replace=TRUE) else xid <- 1:nrow(Z)
if (!is.null(id)) xid <- id
rr <- rr[xid]
Z <- Z[xid,,drop=FALSE]
model <- "fg"
if (cox$prop.odds==TRUE) model <- "logistic"
}# }}}
if (inherits(cox,"phreg"))
{# {{{
p <- length(cox$coef)
if (is.null(Z)) {
Z <- cox$model.frame[,-1,drop=FALSE]
if (cox$nstrata>1) {
ms <- match(cox$strata.name,names(Z))
stratname <- substring(cox$strata.name,8,nchar(cox$strata.name)-1)
strata <- cox$strata
Z <- Z[,-ms,drop=FALSE]
}
}
nz <- ncol(Z)
if (fixZ) Z <- Z else Z <- Z[,names(cox$coef),drop=FALSE]
lrr <- as.matrix(Z) %*% cox$coef
cumhazard <- rbind(c(0,0),cox$cumhaz)
rr <- exp(lrr)
if (drawZ==TRUE) xid <- sample(1:nrow(Z),n,replace=TRUE) else xid <- 1:nrow(Z)
if (!is.null(id)) xid <- id
if (cox$nstrata>1) stratid <- strata[xid] else stratid <- NULL
rr <- rr[xid]
Z <- Z[xid,,drop=FALSE]
if (cox$nstrata>1) {
Z[,stratname] <- stratid
}
model <-c(class(cox),is.null(cox$propodds))
## if (cox$prop.odds==TRUE) cox$model <- "logistic"
}# }}}
if (inherits(cox,"crr"))
{# {{{
p <- length(cox$coef)
if (is.null(Z)) stop("must give covariates for crr simulations\n");
nz <- ncol(Z)
rownames(Z) <- NULL
if (fixZ) Z <- Z else Z <- Z[,names(cox$coef),drop=FALSE]
lrr <- as.matrix(Z) %*% cox$coef
cumhazard <- rbind(c(0,0),cbind(cox$uftime,cumsum(cox$bfitj)))
rr <- exp(lrr)
if (drawZ==TRUE) xid <- sample(1:nrow(Z),n,replace=TRUE) else xid <- 1:nrow(Z)
if (!is.null(id)) xid <- id
rr <- rr[xid]
Z <- Z[xid,,drop=FALSE]
model <- "fg"
}# }}}
if (inherits(cox,"comprisk"))
{# {{{
p <- length(cox$gamma)
formula <- attr(cox, "Formula")
### Z1 <- na.omit(model.matrix(cox,data=data))
p <- nrow(cox$gamma)
if (is.null(Z)) {
Z <- na.omit(get_all_vars(formula,data=data))
nz <- ncol(Z)
Z <- Z[,seq(nz-p+1,nz),drop=FALSE]
}
nz <- ncol(Z)
### if (fixZ) Z <- Z else Z <- Z[,names(cox$gamma),drop=FALSE]
lrr <- as.matrix(Z) %*% cox$gamma
### lrr <- as.matrix(Z) %*% cox$gamma
cumhazard <- cox$cum
rr <- exp(lrr)
if (drawZ==TRUE) xid <- sample(1:nrow(Z),n,replace=TRUE) else xid <- 1:nrow(Z)
if (!is.null(id)) xid <- id
rr <- rr[xid]
Z <- Z[xid,,drop=FALSE]
if (!(cox$model %in% c("fg","logistic2"))) stop("musg be fg or logistic2");
model <- cox$model
}# }}}
out <- list(Z=Z,cum=cumhazard,rr=rr,id=xid,model=model)
if (inherits(cox,"phreg"))
if (cox$nstrata>1) out <- c(out,list(strataid=stratid,strataname=stratname))
return(out)
}# }}}
Try the mets package in your browser
Any scripts or data that you put into this service are public.
mets documentation built on May 29, 2024, 3:51 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions read.fit pre.cifs sim.cifsRestrict sim.cifs sim.cif subdist invsubdist simsubdist sim.cause.cox sim.base sim.cox cause.pchazard.sim rcrisks rcrisk simrchaz addCums lin.approx rchaz
Documented in addCums cause.pchazard.sim invsubdist lin.approx pre.cifs rchaz rcrisk rcrisks read.fit sim.base sim.cause.cox sim.cif sim.cifs sim.cifsRestrict sim.cox simrchaz simsubdist subdist
#' Simulation of Piecewise constant hazard model (Cox).
#'
#' Simulates data from piecwise constant baseline hazard that can also be of
#' Cox type. Censor data at highest value of the break points.
#'
#' For a piecewise linear cumulative hazard the inverse is easy to compute with
#' and delayed entry x we compute
#' \deqn{\Lambda^{-1}(\Lambda(x) + E/RR)},
#' where RR are the relative risks and E is exponential with mean 1.
#' This quantity has survival function
#' \deqn{P(T > t | T>x) = exp(-RR (\Lambda(t) - \Lambda(x)))}.
#'
#' @param cumhazard cumulative hazard, or piece-constant rates for periods defined by first column of input.
#' @param rr relative risk for simulations, alternatively when rr=1 specify n
#' @param n number of simulation if rr not given
#' @param entry delayed entry time for simuations.
#' @param cum.hazard specifies wheter input is cumulative hazard or rates.
#' @param cause name of cause
#' @param extend to extend piecewise constant with constant rate. Default is average rate over time from cumulative (when TRUE), if numeric then uses given rate.
#' @author Thomas Scheike
#' @keywords survival
#' @examples
#'
#' chaz <- c(0,1,1.5,2,2.1)
#' breaks <- c(0,10, 20, 30, 40)
#' cumhaz <- cbind(breaks,chaz)
#' n <- 100
#' X <- rbinom(n,1,0.5)
#' beta <- 0.2
#' rrcox <- exp(X * beta)
#'
#' pctime <- rchaz(cumhaz,n=1000,cum.hazard=FALSE)
#' pctimecox <- rchaz(cumhaz,rrcox,cum.hazard=FALSE)
#'
#' @export
#' @aliases simrchaz addCums lin.approx
rchaz <- function(cumhazard,rr,n=NULL,entry=NULL,cum.hazard=TRUE,cause=1,extend=FALSE)
{# {{{
### cumh=cbind(breaks,rates), first rate is 0 if cumh=FALSE
### cumh=cbind(breaks,cumhazard) if cumh=TRUE
### if (is.null(rr) & is.null(n)) stop("must give rr or n\n");
if (!is.null(n)) rr <- rep(1,n)
n <- length(rr)
breaks <- cumhazard[,1]
rates <- cumhazard[,2][-1]
mm <- tail(breaks,1)
if (cum.hazard==FALSE) {
cumh <- cumsum(c(0,diff(breaks)*rates))
cumhazard <- cbind(breaks,cumh)
} else cumh <- cumhazard[,2]
ttt <- rexp(n)/rr
if (cumhazard[1,2]>0) { ## start cumulative hazard with a 0
### warning("Safest to start with cumulative hazard 0 to avoid problems\n");
cumhazard <- rbind(c(0,0),cumhazard)
cumh <- c(0,cumh)
}
###
if (!is.null(entry)) {
if (length(entry)==1) entry <- rep(entry,n) else entry <- entry
cumentry <- lin.approx(entry,cumhazard,x=1)
} else { entry <- cumentry <- rep(0,n) }
###
ttte <- ttt+cumentry
rrx <- lin.approx(ttte,cumhazard,x=-1)
rrx <- ifelse(rrx>mm,mm,rrx)
status <- rep(0,n)
status <- ifelse(rrx<mm,cause,status)
tcum <- tail(cumhazard,1)
extend.rate <- NULL
if (is.logical(extend))
if (extend) extend.rate <- tcum[2]/tcum[1]
if (is.numeric(extend)) { extend.rate <- extend; extend <- TRUE;}
if (extend) {
whoe <- which(status==0)
nn <- length(whoe)
if (nn>0) {
textend <- rexp(nn)/(rr[whoe]*extend.rate)
rrx[whoe] <- mm+textend
status[whoe] <- 1
}
}
dt <- data.frame(entry=entry,time=rrx,status=status,rr=rr)
colnames(dt) <- c("entry","time","status","rr");
attr(dt,"cumhaz") <- cumhazard
attr(dt,"extend.rate") <- extend.rate
return(dt)
}# }}}
#' @export
lin.approx <- function(x2,xfx,x=1)
{# {{{
### x=1 gives f(x2)
### x=-1 gives f^-1(x2)
breaks <- xfx[,x]
fx <- xfx[,-x]
ri <- fast.approx(breaks,x2,type="left")
maxindex <- which(ri==length(breaks))
rip1 <- ri+1
rip1[maxindex] <- length(breaks)
rrr <- (x2-breaks[ri])/(breaks[rip1]-breaks[ri])
rrr[maxindex] <- 0
res <- rrr*(fx[rip1]-fx[ri])+fx[ri]
res[is.na(res)] <- tail(fx,1)
return(res)
}# }}}
#' @export
addCums <- function(cumB,cumA,max=NULL)
{# {{{
## take times
times <- sort(unique(c(cumB[,1],cumA[,1])))
if (!is.null(max)) times <- times[times<max]
cumBjx <- lin.approx(times,cumB,x=1)
cumAjx <- lin.approx(times,cumA,x=1)
cumBA <- cumBjx+cumAjx
return(cbind(times,cumBA))
}# }}}
#' @export
simrchaz <- function(cumhazard,rr,n=NULL,cens=NULL,rrc=NULL,...)
{# {{{
### adds censoring to to rchaz call
if (!is.null(n)) rr <- rep(1,n)
n <- length(rr)
ptt <- rchaz(cumhazard,rr,...)
if (is.null(rrc)) {
rrc <- rep(1,n)
}
if (is.matrix(cens)) {
pct <- rchaz(cens,rrc,...)
pct <- pct$time
} else {
if (is.numeric(cens)) pct<- rexp(n)/cens else {
chaz <-sum(ptt$status)/sum(ptt$time) ## hazard averate T haz
pct<- rexp(n)/chaz
}
}
dt <- data.frame(time=pmin(ptt$time,pct), status=ifelse(ptt$time<pct,ptt$status,0))
return(dt)
}# }}}
#' Simulation of Piecewise constant hazard models with two causes (Cox).
#'
#' Simulates data from piecwise constant baseline hazard that can also be of
#' Cox type. Censor data at highest value of the break points for either of the
#' cumulatives.
#'
#' @param cumhaz1 cumulative hazard of cause 1
#' @param cumhaz2 cumulative hazard of cause 1
#' @param rr1 number of simulations or vector of relative risk for simuations.
#' @param rr2 number of simulations or vector of relative risk for simuations.
#' @param n number of simulation if rr not given
#' @param cens to censor further , rate or cumumlative hazard
#' @param rrc retlativ risk for censoring.
#' @param ... arguments for rchaz
#' @author Thomas Scheike
#' @keywords survival
#' @examples
#' library(mets); data(bmt); library(survival)
#'
#' cox1 <- phreg(Surv(time,cause==1)~tcell+platelet,data=bmt)
#' cox2 <- phreg(Surv(time,cause==2)~tcell+platelet,data=bmt)
#'
#' X1 <- bmt[,c("tcell","platelet")]
#' n <- 100
#' xid <- sample(1:nrow(X1),n,replace=TRUE)
#' Z1 <- X1[xid,]
#' Z2 <- X1[xid,]
#' rr1 <- exp(as.matrix(Z1) %*% cox1$coef)
#' rr2 <- exp(as.matrix(Z2) %*% cox2$coef)
#'
#' d <- rcrisk(cox1$cum,cox2$cum,rr1,rr2)
#' dd <- cbind(d,Z1)
#'
#' scox1 <- phreg(Surv(time,status==1)~tcell+platelet,data=dd)
#' scox2 <- phreg(Surv(time,status==2)~tcell+platelet,data=dd)
#' par(mfrow=c(1,2))
#' plot(cox1); plot(scox1,add=TRUE)
#' plot(cox2); plot(scox2,add=TRUE)
#' cbind(cox1$coef,scox1$coef,cox2$coef,scox2$coef)
#'
#' @export
#' @aliases cause.pchazard.sim rcrisks
rcrisk <-function(cumhaz1,cumhaz2,rr1,rr2,n=NULL,cens=NULL,rrc=NULL,...)
{#'# {{{
if (!is.null(n)) { rr1 <- rep(1,n); rr2 <- rep(1,n) }
n <- length(rr1);
if (missing(rr2)) rr2 <-rep(1,n)
ptt1 <- rchaz(cumhaz1,rr1,cum.hazard=TRUE,...)
ptt2 <- rchaz(cumhaz2,rr2,cum.hazard=TRUE,...)
ptt <- data.frame(time=pmin(ptt1$time,ptt2$time),status=ifelse(ptt1$time<=ptt2$time,ptt1$status,ptt2$status*2),entry=ptt1$entry)
if (!is.null(cens)) {
if (is.null(rrc)) rrc <- rep(1,n)
if (is.matrix(cens)) {
pct <- rchaz(cens,rrc,...)
pct <- pct$time
} else {
if (is.numeric(cens)) pct<- rexp(n)/cens else {
chaz <-sum(ptt1$status+ptt2$status)/sum(ptt1$time+ptt2$time) ## hazard averate T haz
pct<- rexp(n)/chaz
}
}
ptt <- data.frame(time=pmin(ptt$time,pct),status=ifelse(ptt$time<pct,ptt$status,0),entry=ptt$entry)
}
return(ptt)
}# }}}
#' @export
rcrisks <-function(cumhazs,rrs,n=NULL,cens=NULL,rrc=NULL,entry=NULL,causes=NULL,...)
{#'# {{{
status <- NULL
if (!is.null(n)) rrs <- matrix(1,n,length(cumhazs))
n <- nrow(rrs);
if (!is.null(cens)) cens <- list(cens)
if (!is.null(cens)) {
if (is.null(rrc)) rrc <- rep(1,n);
rrs <- cbind(rrs,rrc)
}
cumss <- c(cumhazs,cens)
for (i in seq_along(cumss)) cumss[[i]] <- rbind(0,cumss[[i]])
cumss <- extendCums(cumss,NULL)
## first time
ptt <- rchaz(cumss[[1]],rrs[,1],entry=entry)
## other times and always min
if (length(cumss)>=2)
for (i in 2:length(cumss)) {
cum <- cumss[[i]]
ptt1 <- rchaz(cumss[[i]],rrs[,i],entry=entry)
ptt <- data.frame(time=pmin(ptt$time,ptt1$time),entry=ptt1$entry,
status=ifelse(ptt$time<=ptt1$time,ptt$status,ptt1$status*i))
}
if (!is.null(cens)) ptt <- dtransform(ptt,status=0,status==length(cumss))
if (!is.null(causes)) {
ptt$status <- c(0,causes)[ptt$status+1]
}
return(ptt)
}# }}}
#' @export
cause.pchazard.sim<-function(cumhaz1,cumhaz2,rr1,rr2,cens=NULL,rrc=NULL,...)
{
rcrisk(cumhaz1,cumhaz2,rr1,rr2,cens=NULL,rrc=NULL,...)
}
#' Simulation of output from Cox model.
#'
#' Simulates data that looks like fit from Cox model. Censor data automatically
#' for highest value of the event times by using cumulative hazard.
#'
#' @param cox output form coxph or cox.aalen model fitting cox model.
#' @param n number of simulations.
#' @param data to extract covariates for simulations (draws from observed
#' covariates).
#' @param cens specifies censoring model, if "is.matrix" then uses cumulative
#' hazard given, if "is.scalar" then uses rate for exponential, and if not
#' given then takes average rate of in simulated data from cox model.
#' @param rrc possible vector of relative risk for cox-type censoring.
#' @param entry delayed entry variable for simulation.
#' @param ... arguments for rchaz, for example entry-time
#' @author Thomas Scheike
#' @keywords survival
#' @examples
#'
#' data(sTRACE)
#' cox <- coxph(Surv(time,status==9)~vf+chf+wmi,data=sTRACE)
#' sim1 <- sim.cox(cox,1000,data=sTRACE)
#' cc <- coxph(Surv(time,status)~vf+chf+wmi,data=sim1)
#' cbind(cox$coef,cc$coef)
#'
#' cor(sim1[,c("vf","chf","wmi")])
#' cor(sTRACE[,c("vf","chf","wmi")])
#'
#' cox <- phreg(Surv(time, status==9)~vf+chf+wmi,data=sTRACE)
#' sim3 <- sim.cox(cox,1000,data=sTRACE)
#' cc <- phreg(Surv(time, status)~vf+chf+wmi,data=sim3)
#' cbind(cox$coef,cc$coef)
#' plot(cox,se=TRUE)
#' plot(cc,add=TRUE,col=2)
#'
#' cox <- phreg(Surv(time,status==9)~strata(chf)+vf+wmi,data=sTRACE)
#' sim3 <- sim.cox(cox,100,data=sTRACE)
#' cc <- phreg(Surv(time, status)~strata(chf)+vf+wmi,data=sim3)
#' cbind(cox$coef,cc$coef)
#' plot(cox)
#' plot(cc,add=TRUE,col=2)
#'
#' @aliases read.fit sim.base simulate.cox
#' @export sim.cox
#' @usage sim.cox(cox,n,data=NULL,cens=NULL,rrc=NULL,entry=NULL,...)
sim.cox <- function(cox,n,data=NULL,cens=NULL,rrc=NULL,entry=NULL,...)
{# {{{
### cumh=cbind(breaks,rates), first rate is 0 if cumh=FALSE
### cumh=cbind(breaks,cumhazard) if cumh=TRUE
des <- read.fit(cox,n,data=data,...)
Z <- des$Z; cumhazard <- des$cum; rr <- des$rr;
if (!inherits(cox,"phreg")) {
if (cumhazard[1,2]>0) cumhazard <- rbind(c(0,0),cumhazard)
ptt <- rchaz(cumhazard,rr,entry=entry)
ptt <- cbind(ptt,Z)
} else {
ptt <- data.frame()
stratj <- cox$strata[cox$jumps]
if (cox$nstrata>1) {
for (j in 0:(cox$nstrata-1)) {
cumhazardj <- rbind(c(0,0),cox$cumhaz[stratj==j,])
if (!is.null(entry)) entry <- entry[des$strataid==j]
pttj <- rchaz(cumhazardj,rr[des$strataid==j],entry=entry)
Zj <- Z[des$strataid==j,,drop=FALSE]
pttj <- cbind(pttj,Zj)
ptt <- rbind(ptt,pttj)
}
} else {
if (cumhazard[1,2]>0) cumhazard <- rbind(c(0,0),cumhazard)
ptt <- rchaz(cumhazard,rr,entry=entry)
ptt <- cbind(ptt,Z)
}
}
if (!is.null(cens)) {
if (is.null(rrc)) rrc <- rep(1,n)
if (is.matrix(cens)) {
pct <- rchaz(cens,rrc,entry=entry)
pct <- pct$time
}
else {
if (is.numeric(cens)) pct<- rexp(n)/cens else {
chaz <-sum(ptt$status)/sum(ptt$time) ## hazard averate T haz
pct<- rexp(n)/chaz
if (!is.null(entry)) pct <- entry + pct
}
}
ptt$time <- pmin(ptt$time,pct)
ptt$status <- ifelse(ptt$time<pct,ptt$status,0)
}
attr(ptt,"id") <- des$id
return(ptt)
}# }}}
#' @export sim.base
#' @usage sim.base(cumhaz,n,data=NULL,strata=NULL,cens=NULL,entry=NULL,...)
sim.base <- function(cumhaz,n,data=NULL,strata=NULL,cens=NULL,entry=NULL,...)
{# {{{
### cumh=cbind(breaks,rates), first rate is 0 if cumh=FALSE
### cumh=cbind(breaks,cumhazard) if cumh=TRUE
if (!is.null(strata)) {
us <- unique(strata)
nus <- length(us)
if (length(n)!=nus) warning("n must be given for each strata as vector \n")
if (length(n)!=nus) n <- rep(round(n/nus),nus)
}
nt <- 0
if (is.null(strata)) {
if (cumhaz[1,2]>0) cumhaz <- rbind(c(0,0),cumhaz)
ptt <- rchaz(cumhaz,n=n,entry=entry)
nt <- n
} else {
ptt <- c()
nstrata <- length(unique(strata))
for (i in seq_along(unique(strata))) {
j <- unique(strata)[i]
cumhazardj <- rbind(c(0,0),cumhaz[strata==j,])
if (!is.null(entry)) entry <- entry[strata==j]
ns <- n[i]
nt <- nt+ns
pttj <- cbind(rchaz(cumhazardj,n=ns,entry=entry),j)
colnames(pttj)[5] <- "strata"
ptt <- rbind(ptt,pttj)
}
}
if (!is.null(cens)) {
if (is.matrix(cens)) {
pct <- rchaz(cens,n=nt,entry=entry)$time
}
else {
pct<- rexp(nt)/(cens)
}
ptt$time <- pmin(ptt$time,pct)
ptt$status <- ifelse(ptt$time<pct,ptt$status,0)
}
return(ptt)
}# }}}
#' Simulation of cause specific from Cox models.
#'
#' Simulates data that looks like fit from cause specific Cox models.
#' Censor data automatically. When censoring is given in the list of causes this
#' will give censoring that looks like the data. Covariates are drawn from data-set
#' with replacement. This gives covariates like the data.
#'
#'
#' @param coxs list of cox models.
#' @param n number of simulations.
#' @param data to extract covariates for simulations (draws from observed
#' covariates).
#' @param cens specifies censoring model, if NULL then only censoring for
#' each cause at end of last event of this type.
#' if "is.matrix" then uses cumulative.
#' hazard given, if "is.scalar" then uses rate for exponential, and if not
#' given then takes average rate of in simulated data from cox model.
#' But censoring can also be given as a cause.
#' @param rrc possible vector of relative risk for cox-type censoring.
#' @param ... arguments for rchaz, for example entry-time
#' @author Thomas Scheike
#' @keywords survival
#' @examples
#' nsim <- 100
#'
#' data(bmt)
#' # coxph
#' cox1 <- coxph(Surv(time,cause==1)~tcell+platelet,data=bmt)
#' cox2 <- coxph(Surv(time,cause==2)~tcell+platelet,data=bmt)
#' coxs <- list(cox1,cox2)
#' dd <- sim.cause.cox(coxs,nsim,data=bmt)
#' scox1 <- coxph(Surv(time,status==1)~tcell+platelet,data=dd)
#' scox2 <- coxph(Surv(time,status==2)~tcell+platelet,data=dd)
#' cbind(cox1$coef,scox1$coef)
#' cbind(cox2$coef,scox2$coef)
#'
#' data(bmt)
#' cox1 <- phreg(Surv(time,cause==1)~tcell+platelet,data=bmt)
#' cox2 <- phreg(Surv(time,cause==2)~tcell+platelet,data=bmt)
#' coxs <- list(cox1,cox2)
#' dd <- sim.cause.cox(coxs,nsim,data=bmt)
#' scox1 <- phreg(Surv(time,status==1)~tcell+platelet,data=dd)
#' scox2 <- phreg(Surv(time,status==2)~tcell+platelet,data=dd)
#' cbind(cox1$coef,scox1$coef)
#' cbind(cox2$coef,scox2$coef)
#' par(mfrow=c(1,2))
#' plot(cox1); plot(scox1,add=TRUE);
#' plot(cox2); plot(scox2,add=TRUE);
#'
#' cox1 <- phreg(Surv(time,cause==1)~strata(tcell)+platelet,data=bmt)
#' cox2 <- phreg(Surv(time,cause==2)~strata(tcell)+platelet,data=bmt)
#' coxs <- list(cox1,cox2)
#' dd <- sim.cause.cox(coxs,nsim,data=bmt)
#' scox1 <- phreg(Surv(time,status==1)~strata(tcell)+platelet,data=dd)
#' scox2 <- phreg(Surv(time,status==2)~strata(tcell)+platelet,data=dd)
#' cbind(cox1$coef,scox1$coef)
#' cbind(cox2$coef,scox2$coef)
#' par(mfrow=c(1,2))
#' plot(cox1); plot(scox1,add=TRUE);
#' plot(cox2); plot(scox2,add=TRUE);
#'
#' @export sim.cause.cox
#' @usage sim.cause.cox(coxs,n,data=NULL,cens=NULL,rrc=NULL,...)
sim.cause.cox <- function(coxs,n,data=NULL,cens=NULL,rrc=NULL,...)
{# {{{
### cumh=cbind(breaks,rates), first rate is 0 if cumh=FALSE
### cumh=cbind(breaks,cumhazard) if cumh=TRUE
if (!is.list(coxs)) stop("Cox models in list form\n");
ptt <- sim.cox(coxs[[1]],n,data=data)
simcovs <- ptt[,(5:ncol(ptt))]
i=2
if (length(coxs)>=2)
for (i in 2:length(coxs)) {
cox <- coxs[[i]]
ptt1 <- sim.cox(coxs[[i]],n,data=data,id=attr(ptt,"id"))
Zi <- ptt1[,(5:ncol(ptt1))]
pm <- match(names(Zi),names(simcovs))
Zi <- Zi[,-pm]
simcovs <- cbind(simcovs,Zi)
dt <- data.frame(time=pmin(ptt$time,ptt1$time),
status=ifelse(ptt$time<=ptt1$time,ptt$status,ptt1$status*i))
}
dt <- cbind(dt,simcovs)
if (!is.null(cens)) {
if (is.matrix(cens)) {
pct <- rchaz(cens,rrc)
pct <- pct$time
}
else {
if (is.numeric(cens)) pct<- rexp(n)/cens else {
chaz <-sum(ptt$status)/sum(ptt$time) ## hazard averate T haz
pct<- rexp(n)/chaz
}
}
dt <- cbind(data.frame(time=pmin(ptt$time,pct),status=ifelse(ptt$time<pct,ptt$status,0)),simcovs)
}
return(dt)
}# }}}
#' @export
simsubdist <- function(cumhazard,rr,n=NULL,entry=NULL,type="cloglog",startcum=c(0,0),...)
{# {{{
## Fine-Gray model cloglog F1= 1-exp(-cum(t)*rr)
## logistic F1= cum(t)*rr/(1+cum(t)*rr)
## identity F1= cum(t)*rr
## rr=exp(X^t beta)
if (!is.null(n)) rr <- rep(1,n)
entry=NULL
logit <- function(p) log(p/(1-p))
breaks <- cumhazard[,1]
rates <- cumhazard[,2][-1]
mm <- tail(breaks,1)
cumh <- cumhazard[,2]
n <- length(rr)
if (cumh[1]>0) {
cumh <- c(startcum[2],cumh)
breaks <- c(startcum[1],breaks)
}
if (type=="cloglog") {
F1tau <- 1-exp(-tail(cumh,1)*rr)
ttt <- -log(1-runif(n)*F1tau)/rr
} else if (type=="logistic") {
F1tau <- tail(cumh,1)*rr/(1+tail(cumh,1)*rr)
v <- runif(n)*F1tau
ttt <- exp(logit(v))/rr;
} else if (type=="identity" | type=="cif") {
F1tau <- tail(cumh,1)
ttt <- runif(n)*F1tau
## rr only affects binomial draw
F1tau <- F1tau*rr
} else stop(" cloglog or logistic or give function (fun=) \n");
###
entry <- cumentry <- rep(0,n)
ttte <- ttt+cumentry
rrx <- lin.approx(ttt,cbind(breaks,cumh),x=-1)
timecause <- rrx
###
rrx <- ifelse(rrx>mm,mm,rrx)
status <- rbinom(n,1,F1tau)
rrx[status==0] <- mm
dt <- data.frame(entry=entry,time=rrx,status=status,rr=rr,F1tau=F1tau,timecause=timecause)
attr(dt,"cumhaz") <- cumhazard
return(dt)
}# }}}
#' @export
invsubdist <- function(F1,u,entry=NULL,cond=1,ptrunc=NULL)
{# {{{
### computes inverse subdist F1^{-1}(u)
### computes inverse subdist F1^{-1}(u)
F1max <- tail(F1[,2],1)
mmax <- tail(F1[,1],1)
if (!is.null(entry)) F1entry <- subdist(F1,entry)[,2];
if (!is.null(entry)) if (is.null(ptrunc)) ptrunc <- 1- F1entry
if (cond==0 & !is.null(entry)) u<-u*ptrunc+F1entry
if (cond==1 & !is.null(entry)) u<-u*(F1max-F1entry)+F1entry
if (cond==1 & is.null(entry)) u<-u*F1max
rrx <- lin.approx(u,F1,x=-1)
status <- rep(1,length(u))
status[u>=F1max] <- 0
dt <- data.frame(time=rrx,status=status,y=u)
if (!is.null(entry)) dt <- cbind(dt,entry)
attr(dt,"F1") <- F1
attr(dt,"Fmax") <- F1max
return(dt)
}# }}}
#' @export
subdist <- function(F1,times)
{# {{{
rrx <- lin.approx(times,F1,x=1)
dt <- cbind(times,rrx)
colnames(dt) <- c("time","subdist")
return(dt)
}# }}}
#' Simulation of output from Cumulative incidence regression model
#'
#' Simulates data that looks like fit from fitted cumulative incidence model
#'
#' @param cif output form prop.odds.subdist or ccr (cmprsk), can also call invsubdist with
#' with cumulative and linear predictor
#' @param n number of simulations.
#' @param data to extract covariates for simulations (draws from observed
#' covariates).
#' @param Z to use these covariates for simulation rather than drawing new ones.
#' @param drawZ to random sample from Z or not
#' @param cens specifies censoring model, if "is.matrix" then uses cumulative
#' hazard given, if "is.scalar" then uses rate for exponential, and if not
#' given then takes average rate of in simulated data from cox model.
#' @param rrc possible vector of relative risk for cox-type censoring.
#' @param cumstart to start cumulatives at time 0 in 0.
#' @param ... arguments for invsubdist
#' @author Thomas Scheike
#' @keywords survival
#' @examples
#' data(bmt)
#'
#' scif <- cifreg(Event(time,cause)~tcell+platelet+age,data=bmt,cause=1,prop=NULL)
#' summary(scif)
#' plot(scif)
#' ################################################################
#' # simulating several causes with specific cumulatives
#' ################################################################
#'
#' cif1 <- cifreg(Event(time,cause)~tcell+age,data=bmt,cause=1,prop=NULL)
#' cif2 <- cifreg(Event(time,cause)~tcell+age,data=bmt,cause=2,prop=NULL)
#' # dd <- sim.cifsRestrict(list(cif1,cif2),200,data=bmt)
#
#' dd <- sim.cifs(list(cif1,cif2),200,data=bmt)
#' scif1 <- cifreg(Event(time,cause)~tcell+age,data=dd,cause=1)
#' scif2 <- cifreg(Event(time,cause)~tcell+age,data=dd,cause=2)
#'
#' par(mfrow=c(1,2))
#' plot(cif1); plot(scif1,add=TRUE,col=2)
#' plot(cif2); plot(scif2,add=TRUE,col=2)
#' @export sim.cif
#' @aliases sim.cif sim.cifs subdist pre.cifs sim.cifsRestrict simsubdist invsubdist
#' @usage sim.cif(cif,n,data=NULL,Z=NULL,drawZ=TRUE,cens=NULL,rrc=NULL,cumstart=c(0,0),...)
sim.cif <- function(cif,n,data=NULL,Z=NULL,drawZ=TRUE,cens=NULL,rrc=NULL,cumstart=c(0,0),...)
{# {{{
## also extracts coefficients and baseline from coxph, cox.aalen, phreg
## and uses them assuming that the model is cloglog unless otherwise
if (!inherits(cif,"defined")) {
des <- read.fit(cif,n,data=data,Z=Z,drawZ=drawZ,...)
Z <- des$Z; cumhaz <- des$cum; rr <- des$rr;
id <- des$id
cumhaz <- rbind(cumstart,cumhaz)
znames <- names(Z);
model <- des$model
} else {
cumhaz <- cif$cumhaz
if (is.null(Z)) stop("must give Z")
rr <- exp(as.matrix(Z) %*% cif$coef)
if (is.data.frame(Z)) znames <- names(Z) else znames <- colnames(Z)
model <- cif$model
id <- 1:nrow(Z)
}
if (!inherits(cif,"phreg")) {
if (model=="logistic2" | model=="logistic") ptt <- simsubdist(cumhaz,rr,type="logistic",...) else ptt <- simsubdist(cumhaz,rr,type="cloglog",...)
ptt <- cbind(ptt,Z)
} else { ### phreg class# {{{
ptt <- data.frame()
if (cif$nstrata>1) {
stratj <- cif$strata[cif$jumps]
for (j in 0:(cif$nstrata-1)) {
cumhazardj <- rbind(c(0,0),cif$cumhaz[stratj==j,])
if (model[3])
pttj <- simsubdist(cumhazardj,rr[des$strataid==j],type="cloglog") else
pttj <- simsubdist(cumhazardj,rr[des$strataid==j],type="logistic")
Zj <- Z[des$strataid==j,,drop=FALSE]
pttj <- cbind(pttj,Zj)
ptt <- rbind(ptt,pttj)
ptt <- rbind(ptt,pttj)
}
} else {
if (model[3]) ptt <- simsubdist(cumhaz,rr,type="cloglog") else
ptt <- simsubdist(cumhaz,rr,type="logistic")
ptt <- cbind(ptt,Z)
}
}# }}}
### adds censoring
if (!is.null(cens)) {# {{{
if (is.null(rrc)) rrc <- rep(1,n)
if (is.matrix(cens)) {
pct <- rchaz(cens,rrc,...)
pct <- pct$time
}
else {
if (is.numeric(cens)) pct<- rexp(n)/cens else {
chaz <-sum(ptt$status)/sum(ptt$time) ## hazard averate T haz
pct<- rexp(n)/chaz
}
}
ptt$time <- pmin(ptt$time,pct)
ptt$status <- ifelse(ptt$time<pct,ptt$status,0)
} # }}}
attr(ptt,"model") <- model
attr(ptt,"id") <- id
attr(ptt,"znames") <- znames
return(ptt)
}# }}}
#' @export sim.cifs
sim.cifs <- function(cifs,n,data=NULL,Z=NULL,cens=NULL,rrc=NULL,max.times=NULL,causes=c(1,2),...)
{# {{{
if (!is.list(cifs)) stop("Cif models in list form\n");
## must consider all models out to last observation times or max.times
cifs <- pre.cifs(cifs,max.times=max.times)
tau <- tail(cifs[[1]]$cum[,1],1)
sim1 <- sim.cif(cifs[[1]],n,data=data,Z=Z)
Z <- sim1[,attr(sim1,"znames")]
if (!inherits(cifs[[2]],"defined")) {
sim2p <- read.fit(cifs[[2]],1,data=data,Z=NULL)
Z2 <- data[attr(sim1,"id"),names(sim2p$Z)]
sim2 <- sim.cif(cifs[[2]],n,data=data,Z=Z2,drawZ=FALSE)
} else {
Z2 <- Z[,attr(cifs[[2]],"znames")]
sim2 <- sim.cif(cifs[[2]],n,data=data,Z=Z2)
}
ptot <- sim1$F1tau+sim2$F1tau
###
rt <- rbinom(n,1,pmin(ptot,1))
rb <- rbinom(n,1,sim1$F1tau/ptot)
cause=ifelse(rb==1,1,2)
time=ifelse(cause==causes[1],sim1$timecause,sim2$timecause)
cause <- rt*cause
time[cause==0] <- tau
ptt <- data.frame(time=time,status=cause,cause=cause,ptot=ptot)
Zcovs <- cbind(Z,Z2)
samecovs <- match(names(Z2),names(Z))
Ze <- Zcovs[,-samecovs]
ptt <- cbind(ptt,Ze)
if (!is.null(cens)) {# {{{
if (is.null(rrc)) rrc <- rep(1,n)
if (is.matrix(cens)) {
pct <- rchaz(cens,rrc,...)
pct <- pct$time
}
else {
if (is.numeric(cens)) pct<- rexp(n)/cens else {
chaz <-sum(ptt$status)/sum(ptt$time) ## hazard averate T haz
pct<- rexp(n)/chaz
}
}
ptt$time <- pmin(ptt$time,pct)
ptt$status <- ifelse(ptt$time<pct,ptt$status,0)
} # }}}
return(ptt)
}# }}}
#' @export sim.cifsRestrict
sim.cifsRestrict <- function(cifs,n,data=NULL,Z=NULL,cens=NULL,rrc=NULL,max.times=NULL,causes=c(1,2),...)
{# {{{
if (!is.list(cifs)) stop("Cif models in list form\n");
## must consider all models out to last observation times or max.times
cifs <- pre.cifs(cifs,max.times=max.times)
tau <- tail(cifs[[1]]$cum[,1],1)
sim1 <- sim.cif(cifs[[1]],n,data=data,Z=Z)
Z <- sim1[,attr(sim1,"znames")]
if (!inherits(cifs[[2]],"defined")) {
sim2p <- read.fit(cifs[[2]],1,data=data,Z=NULL)
Z2 <- data[attr(sim1,"id"),names(sim2p$Z)]
sim2 <- sim.cif(cifs[[2]],n,data=data,Z=Z2,drawZ=FALSE)
} else {
Z2 <- Z[,attr(cifs[[2]],"znames")]
sim2 <- sim.cif(cifs[[2]],n,data=data,Z=Z2)
}
ptot <- sim1$F1tau+sim2$F1tau*(1-sim1$F1tau)
###
rt <- rbinom(n,1,pmin(ptot,1))
rb <- rbinom(n,1,sim1$F1tau/ptot)
cause=ifelse(rb==1,1,2)
time=ifelse(cause==causes[1],sim1$timecause,sim2$timecause)
cause <- rt*cause
time[cause==0] <- tau
ptt <- data.frame(time=time,status=cause,cause=cause,ptot=ptot)
ptt <- cbind(ptt,Z)
samecovs <- match(names(Z2),names(Z))
Ze <- Z2[,-samecovs]
ptt <- cbind(ptt,Ze)
if (!is.null(cens)) {# {{{
if (is.null(rrc)) rrc <- rep(1,n)
if (is.matrix(cens)) {
cum.hazard <- TRUE
pct <- rchaz(cens,rrc,...)
pct <- pct$time
}
else {
if (is.numeric(cens)) pct<- rexp(n)/cens else {
chaz <-sum(ptt$status)/sum(ptt$time) ## hazard averate T haz
pct<- rexp(n)/chaz
}
}
ptt$time <- pmin(ptt$time,pct)
ptt$status <- ifelse(ptt$time<pct,ptt$status,0)
} # }}}
return(ptt)
}# }}}
#' @export pre.cifs
pre.cifs <- function(cifs,n,data=NULL,pprint=FALSE,max.times=NULL,...)
{# {{{
if (!is.list(cifs)) stop("Cif models in list form\n");
## iff needed put cumulative baseline in argument cifs$cum
## extracts maximum time of all cumulatives
maxtimes <- rep(0,length(cifs))
for (i in 1:length(cifs)) {
if (inherits(cifs[[i]],"coxph")) {
stop("Use phreg of mets instead\n");
} else if (inherits(cifs[[i]],"crr")) {
cifs[[i]]$cum <- rbind(c(0,0),cbind(cifs[[i]]$uftime,cumsum(cifs[[i]]$bfitj)))
maxtimes[i] <- max(cifs[[i]]$uftime)
} else if (inherits(cifs[[i]],"phreg")) {
cifs[[i]]$cumhaz <- cifs[[i]]$cumhaz
maxtimes[i] <- max(cifs[[i]]$cumhaz[,1])
} else if (inherits(cifs[[i]],"defined")) {
cifs[[i]]$cumhaz <- cifs[[i]]$cumhaz
maxtimes[i] <- max(cifs[[i]]$cumhaz[,1])
} else {
maxtimes[i] <- max(cifs[[i]]$cum[,1])
}
}
mmtimes <- min(maxtimes)
if (is.null(max.times)) mtimes <- mmtimes else mtimes <- max.times
if (mtimes > mmtimes) {
warning("max.time is further out than max for some cumulatives\n");
cat(" Max times for cumulatives in cifs \n");
print(maxtimes)
}
for (i in 1:length(cifs)) {
cums <- cifs[[i]]$cum
keep <- cums[,1]<mtimes
Fmm <- subdist(cums,mtimes)
cums <- cums[keep,,drop=FALSE]
cums <- rbind(cums,Fmm)
cifs[[i]]$cum <- cums
cifs[[i]]$cumhaz <- cums
if (pprint) {
print(head(cums))
print(tail(cums))
}
}
return(cifs)
}# }}}
## reads a coxph, cox.aalen, phreg, crr, comprisk, prop.odds.subdist object
## and returns cumulative hazard, linear predictor of a resample of size
## n of data, for coxph cox.aalen, comprisk the design Z is constructed
## using the data, but Z can also be given which is needed for crr
## if drawZ is true the covariates in Z are used but multplied coefficient
## based on column names, if fixZ=TRUE the matrix Z is multiplied coefficients
## as is Z %*% coef(x) to form linear predictor
read.fit <- function(cox,n,data=NULL,Z=NULL,drawZ=TRUE,fixZ=FALSE,id=NULL)
{# {{{
if (inherits(cox,"coxph"))
{# {{{
## basehaz(cox,centered=FALSE)[,c(2,1)]
sfit <- survival::survfit(cox, se.fit=FALSE)
zcoef <- ifelse(is.na(coef(cox)), 0, coef(cox))
offset <- sum(cox$means * zcoef)
chaz <- sfit$cumhaz * exp(-offset)
base <- cbind(sfit$time,chaz)
mt <- model.frame(cox)
Y <- model.extract(mt, "response")
if (!inherits(Y, "Surv"))
stop("Response must be a survival object")
if (attr(Y, "type") == "right") {
time <- Y[, "time"];
status <- Y[, "status"]
} else stop("Expected right-censored data.");
if (is.null(Z)) Z <- na.omit(model.matrix(cox))
nn <- colnames(Z)
if (fixZ) Z <- Z else Z <- Z[,names(coef(cox)),drop=FALSE]
rownames(Z) <- NULL
jtime <- sort(time[status==1])
cumhazard <- rbind(c(0,0),Cpred(base,jtime))
rr <- exp( Z %*% matrix(coef(cox),ncol=1))
if (drawZ==TRUE) xid <- sample(1:nrow(Z),n,replace=TRUE) else xid <- 1:nrow(Z)
if (!is.null(id)) xid <- id
rr <- rr[xid]
Z <- Z[xid,,drop=FALSE]
colnames(Z) <- nn
model <- "fg"
}# }}}
if (inherits(cox,"cox.aalen"))
{# {{{
formula <- attr(cox, "Formula")
### Z1 <- na.omit(model.matrix(cox,data=data))
p <- nrow(cox$gamma)
if (is.null(Z)) {
Z <- na.omit(get_all_vars(formula,data=data))
nz <- ncol(Z)
Z <- Z[,seq(nz-p+1,nz)]
}
lrr <- as.matrix(Z) %*% cox$gamma
cumhazard <- cox$cum
rr <- exp(lrr)
if (drawZ==TRUE) xid <- sample(1:nrow(Z),n,replace=TRUE) else xid <- 1:nrow(Z)
if (!is.null(id)) xid <- id
rr <- rr[xid]
Z <- Z[xid,,drop=FALSE]
model <- "fg"
if (cox$prop.odds==TRUE) model <- "logistic"
}# }}}
if (inherits(cox,"phreg"))
{# {{{
p <- length(cox$coef)
if (is.null(Z)) {
Z <- cox$model.frame[,-1,drop=FALSE]
if (cox$nstrata>1) {
ms <- match(cox$strata.name,names(Z))
stratname <- substring(cox$strata.name,8,nchar(cox$strata.name)-1)
strata <- cox$strata
Z <- Z[,-ms,drop=FALSE]
}
}
nz <- ncol(Z)
if (fixZ) Z <- Z else Z <- Z[,names(cox$coef),drop=FALSE]
lrr <- as.matrix(Z) %*% cox$coef
cumhazard <- rbind(c(0,0),cox$cumhaz)
rr <- exp(lrr)
if (drawZ==TRUE) xid <- sample(1:nrow(Z),n,replace=TRUE) else xid <- 1:nrow(Z)
if (!is.null(id)) xid <- id
if (cox$nstrata>1) stratid <- strata[xid] else stratid <- NULL
rr <- rr[xid]
Z <- Z[xid,,drop=FALSE]
if (cox$nstrata>1) {
Z[,stratname] <- stratid
}
model <-c(class(cox),is.null(cox$propodds))
## if (cox$prop.odds==TRUE) cox$model <- "logistic"
}# }}}
if (inherits(cox,"crr"))
{# {{{
p <- length(cox$coef)
if (is.null(Z)) stop("must give covariates for crr simulations\n");
nz <- ncol(Z)
rownames(Z) <- NULL
if (fixZ) Z <- Z else Z <- Z[,names(cox$coef),drop=FALSE]
lrr <- as.matrix(Z) %*% cox$coef
cumhazard <- rbind(c(0,0),cbind(cox$uftime,cumsum(cox$bfitj)))
rr <- exp(lrr)
if (drawZ==TRUE) xid <- sample(1:nrow(Z),n,replace=TRUE) else xid <- 1:nrow(Z)
if (!is.null(id)) xid <- id
rr <- rr[xid]
Z <- Z[xid,,drop=FALSE]
model <- "fg"
}# }}}
if (inherits(cox,"comprisk"))
{# {{{
p <- length(cox$gamma)
formula <- attr(cox, "Formula")
### Z1 <- na.omit(model.matrix(cox,data=data))
p <- nrow(cox$gamma)
if (is.null(Z)) {
Z <- na.omit(get_all_vars(formula,data=data))
nz <- ncol(Z)
Z <- Z[,seq(nz-p+1,nz),drop=FALSE]
}
nz <- ncol(Z)
### if (fixZ) Z <- Z else Z <- Z[,names(cox$gamma),drop=FALSE]
lrr <- as.matrix(Z) %*% cox$gamma
### lrr <- as.matrix(Z) %*% cox$gamma
cumhazard <- cox$cum
rr <- exp(lrr)
if (drawZ==TRUE) xid <- sample(1:nrow(Z),n,replace=TRUE) else xid <- 1:nrow(Z)
if (!is.null(id)) xid <- id
rr <- rr[xid]
Z <- Z[xid,,drop=FALSE]
if (!(cox$model %in% c("fg","logistic2"))) stop("musg be fg or logistic2");
model <- cox$model
}# }}}
out <- list(Z=Z,cum=cumhazard,rr=rr,id=xid,model=model)
if (inherits(cox,"phreg"))
if (cox$nstrata>1) out <- c(out,list(strataid=stratid,strataname=stratname))
return(out)
}# }}}
Try the mets package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.