TABLE OF CONTENTS

simSurvival/generateevents [ Functions ]

NAME

    generateevents --- Generate random event times

FUNCTION

Generate a set of random event times for a sample, given its size, a single covariate, frailties and regression coefficients. Many baseline hazard specifications are supported, see code comments.

SYNOPSIS

629 generateevents <- function(m, Ji, beta, Ui, Zij, type, params)

INPUTS

    m      number of clusters
    Ji     vector of cluster sizes
    beta   a single true regression coefficient
    Ui     length m vector of "true" frailties
    Zij    length sum(Ji) vector of covariates
    type   type of baseline hazard specification
    params parameters for the baseline hazard

OUTPUTS

    Tij    length sum(Ji) vector of event times

SOURCE

632 {
633     Tij <- rep(0, sum(Ji))
634     Uij <- rep(Ui, Ji)
635     #  Weibull baseline hazard, baseline params$lambda0 and shape params$gamweib
636     if(type == "weibull"){
637         if(!all(c("lambda0", "gamweib")%in%names(params))) stop("Parameters lambda0, gamweib, w not specified for type weibull")
638         lambda0 <- params$lambda0
639         gamweib <- params$gamweib
640         for(ind in 1:sum(Ji)){
641             Tij[ind] <- ((-exp(-beta * Zij[ind]) * log(1 - runif(1)) / 
642                 (lambda0 * Uij[ind]))^(1 / gamweib))
643         }
644     }
645     #  stepfunction baseline hazard, with breakpoints at params$breaks and hazards
646     #  at params$haz
647     if(type == "stepfunction"){
648         breaks <- params$breaks
649         haz <- params$haz
650         if(length(haz) != length(breaks) + 1) stop("Step function params: haz should be one longer than breaks")
651         for(ind in 1:sum(Ji)){
652             accept <- FALSE
653             Tijprop <- 0
654             maxhaz <- max(haz) * Uij[ind] * exp(beta * Zij[ind])
655             #  Use accept-reject sampling
656             while(!accept){
657                 Tijprop <- Tijprop -1 / maxhaz * log(runif(1))
658                 thishaz <- haz[findInterval(Tijprop, breaks) + 1] * Uij[ind] * 
659                     exp(beta * Zij[ind])
660                 if(maxhaz * runif(1) < thishaz){
661                     Tij[ind] <- Tijprop
662                     accept <- TRUE
663                 }
664             }
665         }
666     }
667     #  B-spline baseline hazard. params$b is an object returned by bs(),
668     #  and params$w is a set of weights
669     if(type == "bspline"){
670         b <- params$b
671         w <- params$w
672         rbound <- attr(b, "Boundary.knots")[2]
673         survfn <- bs.survfn(b, w, 1000, rbound)
674         if(survfn>.25) warning(paste("Baseline survival over B - spline support is high:", 
675             format(survfn, digits = 3)))
676         for(ind in 1:sum(Ji)){
677             accept <- FALSE
678             Tijprop <- 0
679             # Accept-reject sampling
680             while(!accept){
681                 maxhaz <- max(w) * Uij[ind] * exp(beta * Zij[ind])
682                 Tijprop <- Tijprop -1 / maxhaz * log(runif(1))
683                 thishaz <- predict(b, min(Tijprop, rbound - 1e-5))%*%w * Uij[ind] * 
684                     exp(beta * Zij[ind])
685                 if(maxhaz * runif(1) < thishaz){
686                     Tij[ind] <- Tijprop
687                     accept <- TRUE
688                 }
689             }
690         }
691     }
692     return(Tij)
693 }