TABLE OF CONTENTS

simSurvival/sim.sample [ Functions ]

NAME

    sim.sample --- main simulation function

FUNCTION

Generate a simulated sample of survival data

SYNOPSIS

788 sim.sample <- function(m = 10, Ji = rep(5, 10), params = NULL)

INPUTS

    m      number of clusters
    Ji     cluster size
    params a list with multiple optional components. If any
           of these is not given, defaults are used.
       beta        "true" regression coefficient
       haz.type    type of baseline hazard (see generateevents)
       haz.params  parameters for baseline hazard
       frail.type  type of frailty distribution (see generaterandom)
       frail.params    parameters for frailty
       Z.type      type of covariate
       Z.params    params for covariate
       C.type      type of censoring process
       C.params    params for censoring

OUTPUTS

    agdata     An A-G data frame containing a simulated sample
    Ui         "true" frailties for the sample
    params     parameters used to generate the sample

SOURCE

791 {
792     if(length(Ji) == 1) Ji <- rep(Ji, m)
793     if(length(Ji) != m) stop("Length of Ji does not match m")
794     params.in <- params
795     #   Default parameters
796     params.default <- list(
797         beta = 1,
798         haz.type = "weibull",
799         haz.params = list(lambda0 = 1, gamweib = 1.8),
800         frail.type = "lognormal",
801         frail.params = list(mu = 1, sigma2=.25),
802         Z.type = "normal",
803         Z.params = list(mu = 0, sigma2 = 1),
804         C.type = "weibull",
805         C.params = list(lambda0 = 1, gamweib = 1.8)
806     )
807     params <- params.default
808     if(!is.null(params.in)){
809         for(n in names(params.in)) eval(parse(text = paste("params$", n, " <- params.in$", 
810             n, sep = "")))
811         if(params.in$haz.type == "bspline" & is.null(params.in$haz.params)){
812             tmax = 3;
813             N = 4
814             b <- bs(0, knots = seq(from = 0, to = tmax, length = N + 1)[ - (N + 1)], 
815                 Boundary.knots = c(0, 2), degree = 3)
816             w <- c(.3, .2, .4, .2, .6, .8, 1) * 4
817             params$haz.params <- list(b = b, w = w)
818         }
819     }
820     # Make covariates
821     Zij <- generaterandom(sum(Ji), params$Z.type, params$Z.params)
822     # Make censoring
823     Cij <- generaterandom(sum(Ji), params$C.type, params$C.params)
824     if(!is.null(params$C.max)) Cij[Cij > params$C.max] <- params$C.max
825     # Make frailties
826     Ui <- generaterandom(m, params$frail.type, params$frail.params)
827     # Make event times
828     Tij <- generateevents(m, Ji, params$beta, Ui, Zij, params$haz.type, params$haz.params)
829     # Make event indicators
830     deltaij <- as.numeric(Tij < Cij)
831     # apply censoring
832     Tij[deltaij == 0] <- Cij[deltaij == 0]
833     
834     # make output data frame
835     agdata <- makeagdata(m, Ji, Tij, deltaij, data.frame(Z = Zij))
836     return(list(agdata = agdata, Ui = Ui, params = params))
837 }