R/generate_spm.R

In BDIB: Bayesian Dynamic Information Borrowing for Longitudinal Data with Informative Dropout

#' Generate longitudinal data from a shared-parameter model (SPM)
#' with MNAR monotone dropout (continuous or binary outcomes)
#'
#' @param n Integer. Sample size (number of subjects).
#' @param t_vec Numeric vector of visit times.
#' @param family Character. Either \code{"gaussian"} for a continuous outcome
#'   or \code{"bernoulli"} for a binary outcome.
#' @param alpha0,alpha1 Numerics. Intercept and time slope in the dropout
#'   model (logit link).
#' @param rho Numeric. Coupling between the shared random effect \eqn{b_i}
#'   and the dropout model; \eqn{\rho \neq 0} induces MNAR.
#' @param sigma_b Numeric (>0). Standard deviation of the subject-specific
#'   random intercept \eqn{b_i \sim \mathcal{N}(0, \sigma_b^2)}.
#' @param theta Numeric. Target effect:
#'   \itemize{
#'     \item \code{family="gaussian"}: change at final visit \eqn{_K}:
#'       \eqn{\theta = \beta_1 t_K} (so the slope is \eqn{\beta_1=\theta/t_K}).
#'     \item \code{family="bernoulli"}: marginal response probability at the
#'       final visit in \eqn{(0,1)} (probit link).
#'   }
#' @param beta0 Numeric. For \code{family="bernoulli"} only: probit intercept
#'   in \eqn{\Phi(\beta_0 + \beta_1 t_k + b_i)}. Ignored for \code{"gaussian"}.
#' @param sigma_e Numeric (>0). For \code{family="gaussian"} only: residual SD.
#'
#' @details
#' \strong{Outcome model (shared-parameter).}
#' \itemize{
#'   \item \emph{Gaussian:} \eqn{Y_{ik}\mid b_i \sim \mathcal{N}(\beta_1 t_k + b_i, \sigma_e^2)},
#'         with \eqn{\beta_1 = \theta/t_K}.
#'   \item \emph{Bernoulli (probit):} \eqn{Y_{ik}\mid b_i \sim \mathrm{Bernoulli}(\Phi(\beta_0 + \beta_1 t_k + b_i))},
#'         with \eqn{\beta_1 = \{\sqrt{1+\sigma_b^2}\,\Phi^{-1}(\theta) - \beta_0\}/t_K}.
#' }
#'
#' \strong{Dropout model (MNAR, monotone).}
#' \deqn{\mathrm{logit}\,\Pr(R_{ik}=0 \mid R_{i,k-1}=1, b_i) = \alpha_0 + \alpha_1 t_k + \rho b_i,}
#' and once dropout occurs (\eqn{R_{ik}=0}), all subsequent \eqn{R_{i,k'}=0} for \eqn{k'>k}.
#'
#' @return A list with:
#' \item{Ymat}{\eqn{n \times K} matrix of outcomes (NA for missing).}
#' \item{Rmat}{\eqn{n \times K} matrix of response indicators (1=observed, 0=missing).}
#' \item{dat}{Long-format \code{data.frame} with columns \code{id}, \code{time}, \code{R}, \code{Y}.}
#' @export
#' @examples
#' set.seed(1000)
#' t_vec <- seq(4, 48, 4)
#' # Continuous outcome
#' out_c <- generate_spm(
#'   n = 50, t_vec = t_vec, family = "gaussian",
#'   alpha0 = -4, alpha1 = 0.05, rho = -1, sigma_b = 1,
#'   theta = -3.2, sigma_e = 1
#' )
#' # Binary outcome
#' out_b <- generate_spm(
#'   n = 200, t_vec = t_vec, family = "bernoulli",
#'   alpha0 = -4, alpha1 = 0.05, rho = -1, sigma_b = 1,
#'   theta = 0.4, beta0 = -2
#' )
generate_spm <- function(
    n, t_vec, family = c("gaussian","bernoulli"),
    alpha0, alpha1, rho, sigma_b,
    theta, beta0 = NA_real_, sigma_e = NA_real_
){
  family <- match.arg(family)
  stopifnot(length(n) == 1L, is.numeric(n), n > 0)
  stopifnot(is.numeric(t_vec), length(t_vec) >= 2L, all(diff(t_vec) > 0))
  stopifnot(is.numeric(alpha0), is.numeric(alpha1), is.numeric(rho))
  stopifnot(is.numeric(sigma_b), sigma_b > 0)
  stopifnot(is.numeric(theta), length(theta) == 1L)

  K  <- length(t_vec)
  tK <- t_vec[K]

  # Family-specific slope and RNG for Y_ik
  if (family == "gaussian") {
    stopifnot(is.numeric(sigma_e), sigma_e > 0)
    slope <- theta / tK
  } else { # bernoulli with probit link; theta is Pr(Y=1) at t_K
    if (!(theta > 0 && theta < 1)) stop("For family='bernoulli', theta must be in (0,1).")
    if (!is.numeric(beta0)) stop("For family='bernoulli', beta0 must be provided.")
    slope <- (sqrt(1 + sigma_b^2) * qnorm(theta) - beta0) / tK
  }

  # Allocate storage
  Y <- matrix(NA_real_, n, K)
  R <- matrix(0L,      n, K)

  # Shared random intercept
  b <- rnorm(n, 0, sigma_b)

  # Simulate outcomes with monotone MNAR dropout
  if (family == "gaussian") {
    for (i in seq_len(n)) {
      missing_flag <- FALSE
      for (k in seq_len(K))  {
        if(!missing_flag){
          # Probability that Y(i,k) is missed
          p_miss <- plogis(alpha0 + alpha1 * t_vec[k] + rho * b[i])
          R[i, k] <- 1-rbinom(1, 1, p_miss)          # generate missingness
          # If observed, draw Y from Normal distribution
          if (R[i, k] == 1)
            Y[i, k] <- rnorm(1, slope * t_vec[k] + b[i], sigma_e) else  missing_flag <- TRUE
        } else {
          R[i, k] <- 0L
          Y[i, k] <- NA_real_
        }
      }
    }
  }else{
    for (i in seq_len(n)) {
      dropped <- FALSE
      for (k in seq_len(K)) {
        if (!dropped) {
          p_miss <- plogis(alpha0 + alpha1 * t_vec[k] + rho * b[i])
          R[i, k] <- 1L - rbinom(1L, 1L, p_miss)
          if (R[i, k] == 1L) {
            Y[i, k] <-rbinom(1, 1, pnorm(beta0+slope * t_vec[k] + b[i]))
          } else {
            dropped <- TRUE
          }
        } else {
          R[i, k] <- 0L
          Y[i, k] <- NA_real_
        }
      }
    }
  }




  list(
    Ymat = Y,
    Rmat = R,
    dat  = data.frame(
      id   = rep(seq_len(n), each = K),
      time = rep(t_vec,      n),
      R    = as.vector(t(R)),
      Y    = as.vector(t(Y))
    )
  )
}