TABLE OF CONTENTS

simulation/generaterecurrent [ Functions ]

NAME

    generaterecurrent --- generate simulated recurrent event times

FUNCTION

Generates recurrent event times from a Weibull hazard with baseline lambda_0, shape gamma.

SYNOPSIS

612 generaterecurrent <- function(m, Ji, Zij, Uij, beta1, beta2, lambda0, gamweib, timedep, Cij)
613 {

INPUTS

    m      number of clusters
    Ji     cluster sizes
    Zij    covariates
    Uij    frailties
    beta1  coefficient for covariate 1
    beta2  coefficient for covariate 2 (if applicable)
    labmda0    weibull baseline
    gamweib    weibull shape
    timedep    boolean to indicate time-dependence

SOURCE

616     # Convert the provided covariates into two matrices of time - dependent covariates
617     if(timedep){
618         Zij1 <- Zij$Zij1
619         Zij1times <- Zij$Zij1times
620         Zij2 <- matrix(0, dim(Zij1)[1], 100)  
621         
622         Rij <- rep(0, sum(Ji) * 100); dim(Rij) <- c(sum(Ji), 100)
623         for(ind in 1:sum(Ji))
624         {
625             timescum <- cumsum(Zij1times[ind, ])
626             timescum[min(which(timescum > Cij[ind]))] <- Inf
627             obstimes <- c(timescum[timescum < Cij[ind]], Cij[ind])
628             Rijcum <- 0;k <- 1
629             while(k < 100){
630                 
631                 hazards <- Uij[ind] * lambda0 * gamweib * (obstimes - Rijcum)^(gamweib - 1) *
632                     exp(beta1 * Zij1[ind, 1:length(obstimes)] + beta2 *
633                     Zij2[ind, 1:length(obstimes)])
634                 maxhaz <- max(hazards, na.rm = TRUE)
635                 accept <- FALSE
636                 # Generate by accept-reject
637                 n <- 0;Rijprop <- 0
638                 while(!accept){
639                     genunif <- runif(1)
640                     Rijprop <- Rijprop - 1 / maxhaz * log(genunif)
641                     thistime <- min(which(timescum > (Rijcum + Rijprop)))
642                     if(thistime == 101) thistime <- 100
643                     thishaz <- Uij[ind] * lambda0 * gamweib * (Rijprop)^(gamweib - 1) *
644                         exp(beta1 * Zij1[ind, thistime] + beta2 * Zij2[ind, thistime])
645                     accunif <- runif(1)
646                     if(accunif * maxhaz < thishaz){
647                         Rij[ind, k] <- Rijprop
648                         Rijcum <- Rijcum + Rijprop
649                         if(Rijcum > Cij[ind]) k <- 100
650                         accept <- TRUE
651                     }
652                     n <- n + 1
653                 }
654             k <- k + 1
655             }   
656         }
657         return(Rij)     
658     }else{
659         if(!is.matrix(Zij) && !is.list(Zij)) {
660             Zij1 <- matrix(rep(Zij, 100), length(Zij), 100)
661             Zij2 <- matrix(0, dim(Zij1)[1], 100)
662         }
663         if(is.list(Zij)){
664             Zij1 <- matrix(rep(Zij$Zij1, 100), length(Zij$Zij1), 100)
665             Zij2 <- matrix(rep(Zij$Zij2, 100), length(Zij$Zij2), 100)
666         }    
667         # Initialize output matrix
668         Rij <- rep(0, sum(Ji) * 100)  # Gap data for times between events
669         dim(Rij) <- c(sum(Ji), 100)
670         # Generate interevent gap times by inversion
671         for(k in 1:100)  Rij[, k] <- ((-exp(-beta1 * Zij1[, k] - beta2 * Zij2[, k]) *
672             log(1 - runif(sum(Ji))) / (lambda0 * Uij))^(1 / gamweib))
673         Rij[Rij == 0] <- 1
674         return(Rij)
675     }
676 }