R/AFTsimrecurrent.R

In sa4khan/AFTjmr: Recurrent event analysis and joint modeling using accelerated failure time (AFT) models

#########################################################
# Data Simulation
#########################################################
sim0.ewaft<-function(t0,xx,kappa,gam,rho,beta){
pred<-sum(xx*beta)
u<-rho*t0*exp(-pred)
repeat{
e<-rexp(1)
r1<-log1p(-exp(gam*log1p(-exp(-u^kappa))))-e
r2<-exp(log1p(-exp(r1))/gam)
r3<-(-log1p(-r2))^(1/kappa)
w<-as.numeric(r3*exp(pred)/rho-t0)
if(is.finite(w)) break
}
return(w)
}
########################################################
sim0.waft<-function(t0,xx,kappa,gam=1,rho,beta){
pred<-sum(xx*beta)
u<-rho*t0*exp(-pred)
repeat{
e<-rexp(1)
w<-(u^kappa+e)^(1/kappa)*exp(pred)/rho-t0
if(is.finite(w)) break
}
return(w)
}
########################################################
sim0.llaft<-function(t0,xx,kappa,gam=1,rho,beta){
pred<-sum(xx*beta)
u<-rho*t0*exp(-pred)
repeat{
e<-rexp(1)
w<-(expm1(log1p(u^kappa)+e))^(1/kappa)*exp(pred)/rho-t0
if(is.finite(w)) break
}
return(w)
}
#############################################
sim1.ewaft<-function(t0,xx,end.time,kappa,gam,rho,beta){
t00<-t0
t<-NULL
repeat{
tt0<-t0
t1<-t0+sim0.ewaft(t0=t0,xx=xx,kappa=kappa,gam=gam,rho=rho,beta=beta)
t<-c(t,t1)
t0<-t1
if(t0>end.time || tt0==t1) break
}
status<-rep(1,length(t))
status[which(t>end.time)]<-0
t[t>end.time]<-end.time
if(length(t)==1){
start<-t00
} else{
start<-c(t00,t[1:(length(t)-1)])
}
data.frame(start=start,stop=t,status=status,t(xx))
}
###################################
sim1.waft<-function(t0,xx,end.time,kappa,gam=1,rho,beta){
t00<-t0
t<-NULL
repeat{
tt0<-t0
t1<-t0+sim0.waft(t0=t0,xx=xx,kappa=kappa,gam=gam,rho=rho,beta=beta)
t<-c(t,t1)
t0<-t1
if(t0>end.time || tt0==t1) break
}
status<-rep(1,length(t))
status[which(t>end.time)]<-0
t[t>end.time]<-end.time
if(length(t)==1){
start<-t00
} else{
start<-c(t00,t[1:(length(t)-1)])
}
data.frame(start=start,stop=t,status=status,t(xx))
}
###################################
sim1.llaft<-function(t0,xx,end.time,kappa,gam=1,rho,beta){
t00<-t0
t<-NULL
repeat{
tt0<-t0
t1<-t0+sim0.llaft(t0=t0,xx=xx,kappa=kappa,gam=gam,rho=rho,beta=beta)
t<-c(t,t1)
t0<-t1
if(t0>end.time || tt0==t1) break
}
status<-rep(1,length(t))
status[which(t>end.time)]<-0
t[t>end.time]<-end.time
if(length(t)==1){
start<-t00
} else{
start<-c(t00,t[1:(length(t)-1)])
}
data.frame(start=start,stop=t,status=status,t(xx))
}
######################################
sim.ewaft<-function(origin,end.time,design.mat,beta,kappa,gam,rho){
n<-length(end.time)
k<-0
sim.data<-NULL
repeat{
k<-k+1
dat0<-sim1.ewaft(t0=origin[k],xx=design.mat[k,],end.time=end.time[k],
   kappa=kappa,gam=gam,rho=rho,beta=beta)
sim.data<-rbind(sim.data,cbind(id=k,dat0))
if(k==n) break
}
return(data.frame(sim.data))
}
######################################
sim.waft<-function(origin,end.time,design.mat,beta,kappa,gam,rho){
n<-length(end.time)
k<-0
sim.data<-NULL
repeat{
k<-k+1
dat0<-sim1.waft(t0=origin[k],xx=design.mat[k,],end.time=end.time[k],
   kappa=kappa,gam=gam,rho=rho,beta=beta)
sim.data<-rbind(sim.data,cbind(id=k,dat0))
if(k==n) break
}
return(data.frame(sim.data))
}
######################################
sim.llaft<-function(origin,end.time,design.mat,beta,kappa,gam,rho){
n<-length(end.time)
k<-0
sim.data<-NULL
repeat{
k<-k+1
dat0<-sim1.llaft(t0=origin[k],xx=design.mat[k,],end.time=end.time[k],
   kappa=kappa,gam=gam,rho=rho,beta=beta)
sim.data<-rbind(sim.data,cbind(id=k,dat0))
if(k==n) break
}
return(data.frame(sim.data))
}
##########################################
#' Simulate recurrent event data
#' @description Simulate recurrent event data for fixed covariates
#' @keywords AFT model, Recurrent event data, Simulation
#' @param surv.formula a one-sided formula of the form \code{~ z1 + z2 + ... + zp}, where 
#'       \code{z1, z2, ..., zp} are baseline covariates for the time-to-event process. 
#' @param surv.par a list of the form \code{list(beta = a vector, logrho = a scalar, 
#'      logkappa = a scalar)} for the time-to-event model parameters, 
#'      where \code{beta} is the vector of regression coefficients for the baseline covariates as specified in
#'      \code{surv.formula} (without an intercept), 
#'      \code{logrho} is the log of the rate parameter
#'      of the time-to-event distribution (Weibull or log-logistic), and
#'      \code{logkappa} is the log of the shape parameter of the time-to-event distribution 
#'      (Weibull or log-logistic). Note that \code{surv.par}
#'      must be a named list as above. See \code{survreg.fit} for parameterizations of the time-to-event distributions. 
#' @param surv.model the AFT model to be used to describe the recurrent event process. Available options are
#'       \code{"weibull"} and \code{"llogistic"}
#'       for Weibull and log-logistic distributions, respectively.
#' @param Data a data frame containing the covariates \code{z1, z2, ..., zp}  for the time-to-event model.
#'        For \code{n} subjects, it must must have \code{n} rows.
#' @param origin the time origin (an \eqn{n} x 1 vector for \eqn{n} subjects; typically, a vector of zeros).
#' @param end.time the end of follow-up time (an \eqn{n} x 1 vector for \eqn{n} subjects).
#' @details Simulate recurrent event data. See Section 2.4 of Cook and Lawless.
#' @return a data frame containing \code{id}, \code{start}, \code{stop}, \code{status} and covariates.
#' @references Cook RJ and Lawless J, The statistical analysis of recurrent events, Springer, 2007.
#' @author Shahedul Khan <khan@math.usask.ca>
#' @seealso \code{\link{survreg.aft}}
#' @examples
#'   n<-100 # number of subjects
#'   Data<-data.frame(z1=rnorm(n),z2=rbinom(n,1,0.5))
#'   origin<-rep(0,n)
#'   end.time<-runif(n,3,3.5)
#'   surv.par<-list(beta=c(0.5,-0.5),logrho=log(0.25),logkappa=log(2))
#'   sim.ll<-sim.rec(surv.formula=~z1+z2,surv.par=surv.par,surv.model="llogistic",
#'       Data=Data,origin=origin,end.time=end.time)
#'   survreg.aft(Formula=c(start,stop,status)~z1+z2,Data=sim.ll,model="llogistic")
#'   sim.w<-sim.rec(surv.formula=~z1+z2,surv.par=surv.par,surv.model="weibull",
#'       Data=Data,origin=origin,end.time=end.time)
#'   survreg.aft(Formula=c(start,stop,status)~z1+z2,Data=sim.w,model="weibull")
#' @export

sim.rec<-function(surv.formula,surv.par,surv.model,Data,origin,end.time){
if(isFALSE(surv.model=="weibull" || surv.model=="llogistic")){
stop("surv.model must be 'weibull' or 'llogistic'")
}
survvar<-all.vars(surv.formula)
if(isFALSE(all(survvar %in% colnames(Data)))){
stop("variable names in 'surv.formula' did not match
     with the variable names in 'Data'")
}
if(length(origin)!=nrow(Data) || length(end.time)!=nrow(Data)){
stop("must be length(origin) = length(end.time) = nrow(Data)")
}
if(any(end.time<=origin)){
stop("end.time must be > origin")
}
surv.par.names<-names(surv.par)
#if(surv.model=="weibull" || surv.model=="llogistic"){
survpar<-c("beta","logrho","logkappa")
if(isFALSE(all(survpar %in% surv.par.names))){
stop("use names to define 'surv.par' list; at least one
  element of 'surv.par' is missing or not recognized")
}
#}
#if(surv.model=="eweibull"){
#survpar<-c("beta","logrho","logkappa","loggamma")
#if(isFALSE(all(survpar %in% surv.par.names))){
#stop("use names to define 'surv.par' list; at least one
#  element of 'surv.par' is missing or not recognized")
#}
#}  
if(length(surv.par$beta)!=length(survvar)){
stop("length of beta in 'surv.par' did not match with 'surv.formula'")
}
design.mat<-matrix(model.matrix(surv.formula,data=Data)[,-1],ncol=length(survvar))
colnames(design.mat)<-survvar
beta<-surv.par$beta
rho<-exp(surv.par$logrho)
kappa<-exp(surv.par$logkappa)
#if(surv.model=="eweibull"){
#gamma<-exp(surv.par$loggamma)
#} else{
#gamma<-NULL
#}
if(surv.model=="weibull") {sim<-sim.waft(origin=origin,end.time=end.time,
   design.mat=design.mat,beta=beta,kappa=kappa,gam=gamma,rho=rho)}
if(surv.model=="llogistic") {sim<-sim.llaft(origin=origin,end.time=end.time,
   design.mat=design.mat,beta=beta,kappa=kappa,gam=gamma,rho=rho)}
#if(model=="eweibull") {sim<-sim.ewaft(origin=origin,end.time=end.time,
#   design.mat=design.mat,beta=beta,kappa=kappa,gam=gamma,rho=rho)}
return(sim)
}