R/generate_data.R

In XuemeiDing/survival19: Deal with Stratified Cox Proportional Hazard Model

#'Generate Data
#'
#'The function generate_data generates example survival data
#'
#'@param true_beta true parameters
#'@param F_pre a unique list of facilities
#'@param gamma true facility effects of each facility, notice length of gamma must be the same as length of F_pre
#'@param mean_sample_size expected sample size of each facility (stratum), a scalar same in all facilities
#'@param Sigma_z1 Var-Cov matrix of design matrix
#'@param censoring_up_bound maximum value that the censoring time can take
#'@param censoring_trunc maximum observation time
#'@param order whether the output data should be ordered by time
#'
#'@return a list of generated data:
#'\item{delta}{event indicator}
#'\item{z}{design matrix}
#'\item{facility}{which facility each patient belongs to}
#'\item{num_facility}{number of facilities (strata)}
#'\item{time}{obseved time for each patient}
#'\item{true_beta}{true parameters}
#'\item{N}{number of patients}
#'
#'@examples
#'surv_dat = generate_data()
#'
#'@export
#'

generate_data = function(true_beta = c(c(0.5, -0.5, 1, -1, 1.5), rep(0, 5)), F_pre = 1:10, gamma = rep(1, 10),
                         mean_sample_size = 50, Sigma_z1 = diag(length(true_beta)), censoring_up_bound = 30,
                         censoring_trunc = 3, order = T){

  p = length(true_beta)
  num_facility = length(F_pre)
  n_f = stats::rpois(num_facility, lambda = mean_sample_size)  #sample size for each facility
  N = sum(n_f)

  facility = rep(F_pre, n_f)
  gamma_subject=rep(gamma, n_f)

  ############ generate data ######################

  z = mvtnorm::rmvnorm(N, mean = rep(0, p), sigma = Sigma_z1)

  U = stats::runif(N, 0, 1)
  pre_time = -log(U)/(gamma_subject * exp(z %*% true_beta))

  pre_censoring = stats::runif(N,1,censoring_up_bound)
  pre_censoring = pre_censoring * (pre_censoring<censoring_trunc) + censoring_trunc * (pre_censoring >= censoring_trunc)

  delta = (pre_censoring >= pre_time) # event indicator

  time = pre_time * (delta == 1) + pre_censoring * (delta == 0)

  ############ order data ######################
  if (order){
    od = order(time)
    delta = delta[od]
    z = z[od, ]
    facility = facility[od]
    time = time[od]
  }

  return(list(delta = delta, z = z, facility = facility, num_facility = num_facility, time = time,
              true_beta = true_beta, N = N))
}