R/generate_data.R

In transmdl: Semiparametric Transformation Models

#' @title  Generate data for simulation.
#'
#' @description   Generate observed event times, covariates and other data used
#'   for simulations in the paper.
#'
#' @return a list containing
#' \tabular{lccl}{
#'   \code{X} \tab\tab\tab design matrix \cr
#'   \code{beta_X} \tab\tab\tab \eqn{X\cdot\beta^T} \cr
#'   \code{Y} \tab\tab\tab observed event time \cr
#' }
#'
#'
#' @examples
#'  generate_data(200,0.5,1,c(0.5,-1))
#'
#'
#' @references Abramowitz, M., and Stegun, I.A. (1972). Handbook of Mathematical
#'   Functions (9th ed.). Dover Publications, New York.
#'+-
#'   Evans, M. and Swartz, T. (2000). Approximating Integrals via Monte Carlo
#'   and Deterministic Methods. Oxford University Press.
#'
#'   Liu, Q. and Pierce, D.A. (1994). A note on Gauss-Hermite quadrature.
#'   Biometrika 81: 624-629.
#'
#'
#' @param n number of subjects
#' @param alpha parameter in transformation function
#' @param rho parameter in baseline cumulative hazard function \eqn{\Lambda(t) =
#'   \rho log(1+t)} assumed in simulation
#' @param beta_true parameter \eqn{\beta}
#' @param now_repeat number of duplication of simulation
#'
#' @details The survival function for \eqn{t} of the \eqn{i}th observation takes
#' the form \deqn{S_{i}(t| X_{i}) = \exp\left\{-H \{\Lambda(t) \exp ( \beta^T
#' X_{i} ) \}\right\}.} The failure time \eqn{T_i} can be generated by \deqn{
#' T_i = \left\{\begin{array}{l l} \exp\{
#' \frac{U^{-\alpha}-1}{\alpha\rho\exp\{\beta^TX_i \}} \}-1& \alpha > 0, \\
#' \exp\{  \frac{-log(U)}{\rho\exp\{\beta^TX_i \}} \}-1, & \alpha = 0.
#' \end{array}\right\} }
#'
#' @export
generate_data = function(n, alpha, rho, beta_true, now_repeat=1)
{

  set.seed(7+700*(now_repeat+1))
  X1 = stats::rbinom(n,1,0.5)
  epsil = stats::rnorm(n)
  X2 = X1 + epsil * (abs(epsil) < 1 ) + (abs(epsil) >= 1 )
  X = matrix(c(X1,X2) , ncol = 2, byrow = FALSE)
  U = stats::runif(n)
  beta_X = c(X %*% beta_true)

  if(alpha == 0)
  {
    t = (exp(-log(U)/(rho * exp(beta_X))) - 1)
  }else
  {
    t=(exp((U^(-alpha) - 1)/(alpha * rho * exp(beta_X))) - 1)
  }
  C = stats::runif(n,2,6)
  Y = pmin(C,t)
  delta_i = ifelse( C >= t, 1, 0)
  return(list(X = X, delta = delta_i, beta_X = beta_X, Y = Y))
}