R/simPregProp.R

In simPreg: Simulate Pregnancy Data with Time-Varying Exposure

#' Simulate distributions of pregnancy outcomes and gestational duration
#'
#' Calculates proportions of pregnancies experiencing spontaneous live birth,
#' non-spontaneous live birth, or late miscarriage/stillbirth by gestational
#' age and timing of exposure, given daily hazard rates and hazard ratios.
#'
#' @param haz.spont.livebirth Numeric vector of daily hazards of spontaneous
#'   live birth. Defaults to the package-provided
#'   `haz.spont.livebirth.default` vector of length 301.
#' @param haz.nonspont.livebirth Numeric vector of daily hazards of
#'   non-spontaneous live birth. Defaults to the package-provided
#'   `haz.nonspont.livebirth.default` vector of length 301.
#' @param haz.late.miscarriage.stillbirth Numeric vector of daily hazards of
#'   late miscarriage/stillbirth. Defaults to the package-provided
#'   `haz.late.miscarriage.stillbirth.default` vector of length 301.
#' @param haz.exposure Numeric vector of daily hazards of exposure. Defaults
#'   to the package-provided `haz.exposure.default` vector of length 301.
#' @param hr.spont.livebirth Numeric vector of hazard ratios for spontaneous
#'   live birth following exposure. Must have the same length as the hazard
#'   vectors. Defaults to `NULL`, corresponding to no effect of exposure on the
#'   hazard.
#' @param hr.nonspont.livebirth Numeric vector of hazard ratios for
#'   non-spontaneous live birth following exposure. Must have the same length
#'   as the hazard vectors. Defaults to `NULL`, corresponding to no effect of
#'   exposure on the hazard.
#' @param hr.late.miscarriage.stillbirth Numeric vector of hazard ratios for
#'   late miscarriage/stillbirth following exposure. Must have the same length
#'   as the hazard vectors. Defaults to `NULL`, corresponding to no effect of
#'   exposure on the hazard.
#'
#' @return A data frame with class `preg.prop` containing the following
#'   columns:
#'   \describe{
#'     \item{GA}{Gestational age in days.}
#'     \item{ExpGA}{Gestational age at exposure. `NA` indicates no exposure.}
#'     \item{Outcome}{Pregnancy outcome: `spontaneous_live_birth`,
#'     `nonspontaneous_live_birth`, or `late_miscarriage_stillbirth`.}
#'     \item{Prop}{Proportion of pregnancies corresponding to the
#'       gestational age, exposure timing, and outcome combination.}
#'   }
#'
#' @examples
#' # Use the package-provided default hazard vectors
#' data.prop <- simPregProp()
#' head(data.prop)
#'
#' # Illustrate an exposure associated with increased risk of preterm live birth
#' hr <- rep(1, 301)
#' hr[1:258] <- 2
#' data.prop <- simPregProp(
#'   hr.spont.livebirth = hr,
#'   hr.nonspont.livebirth = hr
#' )
#'
#' @export
simPregProp <- function(haz.spont.livebirth = haz.spont.livebirth.default,
                        haz.nonspont.livebirth = haz.nonspont.livebirth.default,
                        haz.late.miscarriage.stillbirth = haz.late.miscarriage.stillbirth.default,
                        haz.exposure = haz.exposure.default,
                        hr.spont.livebirth = NULL,
                        hr.nonspont.livebirth = NULL,
                        hr.late.miscarriage.stillbirth = NULL){

  # Collect input arguments
  args <- as.list(environment())

  # Check gestational-age vector length from hazard inputs
  args.haz <- args[startsWith(names(args), "haz.")]
  ga.max <- unique(lengths(args.haz))
  if (length(ga.max) != 1) {
    stop("All hazard vectors must have the same length.")
  }
  if (ga.max < 2) {
    stop("Hazard vectors must contain at least two values.")
  }

  # Check HR vector lengths and set unspecified vectors to 1
  args.hr <- args[startsWith(names(args), "hr.")]
  for (arg.name in names(args.hr)) {
    arg <- args.hr[[arg.name]]
    if (is.null(arg)) {
      assign(arg.name, rep(1, ga.max), envir = environment())
    } else if (length(arg) != ga.max) {
      stop(arg.name, " must have the same length as the hazard vectors.")
    }
  }

  # Update argument list
  args <- as.list(environment())
  args <- args[(startsWith(names(args), "haz.")) | (startsWith(names(args), "hr."))]

  # Check input argument values
  for (arg.name in names(args)) {
    arg <- args[[arg.name]]

    # Check numeric and missing values
    if (!is.numeric(arg)) {
      stop(arg.name, " must be numeric.")
    }
    if (anyNA(arg)) {
      stop(arg.name, " contains missing values.")
    }

    # Check value ranges
    if (startsWith(arg.name, "haz.") &&
        any(arg < 0 | arg > 1)) {
      stop(arg.name, " must be between 0 and 1.")
    }
    if (startsWith(arg.name, "hr.") &&
        any(arg <= 0)) {
      stop(arg.name, " must be > 0.")
    }
  }

  # Check products of hazard rates and hazard ratios
  if (any(haz.spont.livebirth * hr.spont.livebirth > 1)) {
    stop("haz.spont.livebirth * hr.spont.livebirth must be <= 1; one or more values of hr.spont.livebirth must be reduced.")
  }
  if (any(haz.nonspont.livebirth * hr.nonspont.livebirth > 1)) {
    stop("haz.nonspont.livebirth * hr.nonspont.livebirth must be <= 1; one or more values of hr.nonspont.livebirth must be reduced.")
  }
  if (any(haz.late.miscarriage.stillbirth * hr.late.miscarriage.stillbirth > 1)) {
    stop("haz.late.miscarriage.stillbirth * hr.late.miscarriage.stillbirth must be <= 1; one or more values of hr.late.miscarriage.stillbirth must be reduced.")
  }

  # Abbreviate inputs
  haz.slb <- haz.spont.livebirth
  haz.nslb <- haz.nonspont.livebirth
  haz.sb <- haz.late.miscarriage.stillbirth
  haz.exp <- haz.exposure
  hr.slb <- hr.spont.livebirth
  hr.nslb <- hr.nonspont.livebirth
  hr.sb <- hr.late.miscarriage.stillbirth

  # Initialize matrices with rows representing gestational age and columns representing exposure timing
  mat.preg <- matrix(NA, nrow = ga.max, ncol = ga.max)
  mat.sb <- matrix(NA, nrow = ga.max, ncol = ga.max)
  mat.slb <- matrix(NA, nrow = ga.max, ncol = ga.max)
  mat.nslb <- matrix(NA, nrow = ga.max, ncol = ga.max)

  # Initialize first cells
  mat.preg[1,1] <- 1
  mat.sb[1,1] <- mat.preg[1,1]*haz.sb[1]
  mat.slb[1,1] <- (mat.preg[1,1] - mat.sb[1,1])*haz.slb[1]
  mat.nslb[1,1] <- (mat.preg[1,1] - mat.sb[1,1] - mat.slb[1,1])*haz.nslb[1]
  mat.preg[2,2] <- (mat.preg[1,1] - mat.sb[1,1] - mat.slb[1,1] - mat.nslb[1,1])*haz.exp[1]
  mat.preg[2,1] <- mat.preg[1,1] - mat.sb[1,1] - mat.slb[1,1] - mat.nslb[1,1] - mat.preg[2,2]

  # Fill remaining cells
  for (i in 2:nrow(mat.preg)) {
    # Unexposed
    mat.sb[i,1] <- mat.preg[i,1]*haz.sb[i]
    mat.slb[i,1] <- (mat.preg[i,1] - mat.sb[i,1])*haz.slb[i]
    mat.nslb[i,1] <- (mat.preg[i,1] - mat.sb[i,1] - mat.slb[i,1])*haz.nslb[i]
    if (i < nrow(mat.preg)) {
      mat.preg[(i+1),(i+1)] <- (mat.preg[i,1] - mat.sb[i,1] - mat.slb[i,1] - mat.nslb[i,1])*haz.exp[i]
      mat.preg[(i+1),1] <- mat.preg[i,1] - mat.sb[i,1] - mat.slb[i,1] - mat.nslb[i,1] - mat.preg[(i+1),(i+1)]
    }

    # Exposed
    for (j in 2:i) {
      mat.sb[i,j] <- mat.preg[i,j]*haz.sb[i]*hr.sb[i]
      mat.slb[i,j] <- (mat.preg[i,j] - mat.sb[i,j])*haz.slb[i]*hr.slb[i]
      mat.nslb[i,j] <- (mat.preg[i,j] - mat.sb[i,j] - mat.slb[i,j])*haz.nslb[i]*hr.nslb[i]
      if (i < nrow(mat.preg)) {
        mat.preg[(i+1),j] <- mat.preg[i,j] - mat.sb[i,j] - mat.slb[i,j] - mat.nslb[i,j]
      }
    }
  }

  # Store matrices in list
  data.l <- list(spontaneous_live_birth = mat.slb,
                 nonspontaneous_live_birth = mat.nslb,
                 late_miscarriage_stillbirth = mat.sb)

  # Initialize df list of same length
  df.l <- vector("list", length(data.l))

  # Loop through matrices
  for (i in seq_along(data.l)) {

    # Extract matrix
    mat <- data.l[[i]]

    # Get row and column indices of non-zero cells
    indices <- which(!is.na(mat) & mat != 0, arr.ind = TRUE)

    # Create data frame and save in df list
    if (nrow(indices) > 0) {
      df.l[[i]] <- data.frame(GA = indices[,1],
                              ExpGA = indices[, 2] - 1,
                              Outcome = names(data.l)[i],
                              Prop = mat[indices],
                              row.names = NULL)
    }
  }

  # Combine data frames
  df <- do.call(rbind, df.l[!sapply(df.l, is.null)])

  # Set unexposed exposure timing to NA
  df$ExpGA[which(df$ExpGA == 0)] <- NA

  # Remove row names and define class
  rownames(df) <- NULL
  class(df) <- c("preg.prop", class(df))

  # Return
  return(df)
}