Simulation -Discrete Time Hazard model

Question

How we can simulate the survival times for fitting discrete time hazard model?

Answer 1

Assuming you mean a discrete-time logit hazard, complementary log-log hazard, or probit hazard model, you can do this:

1. Establish your model assumptions

Specify a model incorporating at least two categories of assumption:

1. How does $h_{t}$, the discrete-time hazard function depend on time? For example:
   • Is $h_{t}$ constant across time? (I.e. independent of time.)
   • Is $h_{t}$ a linear function of time? (I.e. increases/decreases as a constant rate across time)
   • Is $h_{t}$ a fully discrete function of time (I.e. each period $t$ has it's own conditionally independent intercept term.)
   • Is $h_{t}$ some more complex function of time (E.g., does $h_{t}$ operate as a quadratic function of time? Does $h_{t}$ depend linearly on time, but with one or more discrete effects as specific time periods?)
2. Is $h_{t}$ conditioned on other variables besides time? For example:
   • Is $h_{t}$ unconditioned on non-time predictors? (I.e., is this a baseline model of $h_{t}$?)
   • Is $h_{t}$ conditioned linearly on any predictors?
   • Is $h_{t}$ conditioned non-linearly on any predictors? (E.g., two predictors plus their multiplicative interaction effect, etc.)

2. Prepare for simulation (R sample code, assuming logit hazard, but stay tuned)

1. Build some useful functions

   # create the logit function, and its inverse (over values of x 
   # from 0–1), the logistic function 
   logistic <- function(x) { return(1/(1+exp(-1*x))) }
   logit    <- function(x) { return(log(x/(1-x))) }
   

2. Specify $N$, the total number of individuals (units) at-risk at $T=0$

   N <- 1000   # You need to decide what N is, 1000 is just an example
   

3. Specify $T$, the maximum value of $t$ in your simulated study.

   T <- 20   # You need to decide what T is, 20 is just an example
   

4. Create a data frame called Data starting with the ID variable

   ID <- 1:N
   Data <- data.frame(ID)
   
   # Create conditioning variables (e.g., two dichotomous vars here, but
   # these could be whatever number and kind), and bind these to your 
   # data set
   var1 <- c(rep(0,(N/2)),rep(1,N/2))
   var2 <- c(rep(0,(N/4)),rep(1,(N/2)),rep(0,(N/4)))
   Data <- cbind(Data,var1,var2)
   

5. Format Data as a person-period structure

   expand <- function(x,t) {
     data <- x[rep(1:length(x[,1]), each = t), ]
     rownames(data) <- NULL
     return(data)
     }
   
    Data <- expand(Data, T)
    Data$period <- rep(1:T,N)
        for (t in 1:T) {
          varname <- paste("t",t,sep="")
          Data[,varname] <- as.integer(0 + Data$period == t)
      }
   

6. Prepare effect of time on discrete-time hazard function

   Based on your above model assumptions about how $h_{t}$ relates to time:

   i. Constant baseline hazard:

   # for a nominal hazard of 0.1:
   conslogit   <- logit(0.1)   # logit(.1)   = -2.1972246
   

   ii. Baseline hazard as a linear function of time

   # for a linear effect of time on hazard of 0.05 per 1-unit increase in t:
   linearlogit   <- logit(logistic(0) + 0.05)
   

   iii. Baseline hazard as a fully discrete function of time

   # for a fully discrete time effect on hazard (no. of discrete values
   # should equal T; here there are 20):
   discretehazards <- c(0.11, 0.12, 0.06, 0.20, 0.07, 0.17, 0.05, 0.04, 0.18, 0.06, 0.15, 0.09, 0.16, 0.25, 0.03, 0.15, 0.05, 0.11, 0.08, 0.14)
   discretelogit   <- logit(discretehazards)
   

7. Prepare effects of conditioning variables on $h_{t}$

   # for a nominal effect of var1 of 0.06 increase in hazard:
   var1logit   <- logit(logistic(0) + 0.06)
   #
   # for a nominal hazard effect of var2 of 0.04:
   var2logit   <- logit(logistic(0) + 0.04)
   

3. Simulate your discrete-time survival data

# Example for baseline fully discrete effect of time:
for (t in 1:T) {
  Data$hlogit[Data$period==t] <- logistic(discretelogit[t])
  }

 

# Example for conditional logit hazard model with constant effect of time:
# (The [Data$period==t] on the conditioning variables are probably only
# necessary if your conditioning variables are *time varying*: if the 
# values of var1 or var2 were to change over time *within individuals*.)
for (t in 1:T) {
  Data$hlogit[Data$period==t] <- logistic(conslogit + var1logit*Data$var1[Data$period==t] + var2logit*Data$var2[Data$period==t])
  }

etc.

Data now holds your simulated discrete-time data set!

BONUS: 4. Wait a minute! I am using a discrete-time probit hazard (or complementary log-log) model!

More or less exactly as above, except:

1. For probit hazard

   # The inverse probit function (cumulative distribution function of z)
   invprobit <- function(x) { return(pnorm(x)) }
   probit    <- function(x) { return(qnorm(x)) }
   

2. For complementary log-log hazard

   # The inverse complementary log log function 
   invcloglog <- function(x) { return(1-(exp(-1*exp(x)))) }
   cloglog    <- function(x) { return(log(-1*log(1-x))) }
   

3. Where you have specified effects of time and/or conditioning variable effects using, for example, var1logit <- logit(logistic(0) + 0.06) instead use:

   • Probit: var1probit <- probit(invprobit(0) + 0.06)
   • Complementary log-log: var1cloglog <- cloglog(invcloglog(0) + 0.06)

   and then adjust the final simulation loop to use these variables as appropriate.

PS: The code snippets are based on code from my own simulations. I selectively adapted to answer your question, but there may be errors in my answer, if there are, reach out to me and I will happily assist.

Answer 2

It’s not too difficult. To simulate one survival time:

1. Assume the unit is alive at the start of period/interval 1.

2. Let $h_1$ be the hazard rate for period 1.

3. Take a random draw $x$ from the uniform distribution on the unit interval.

4. If $x \leq h_1$, the unit died in that period/interval. If not, the unit is alive at the start of the next period.

5. If still alive, repeat steps 2-4 the next period, and each subsequent period, until you either simulate a death or you reach some maximum period (after which you can only say that death occurs after that). The max period will be limited by memory and processing time.

Repeat for each unit. Depending on your specific case, you may want to repeat many times for the same unit, repeat across different units (which may have different hazard rates), or both.