R/simulation_functions.R

In gsearly: Creates Group Sequential Trial Designs when Early Outcomes are Available

#-------------------------------------------------------------------------------
# Simulation functions for gsearly
#-------------------------------------------------------------------------------
# 1st July 2026
#-------------------------------------------------------------------------------

#-------------------------------------------------------------------------------
# 5.  .estParams
#-------------------------------------------------------------------------------
.estParams <- function(sim, rdata, odata, s = s, vcovmat = vcovmat, allint) {
  ## Get numbers for simulation number sim
  pick0 <- function(x, rdata, sim) {
    return(rdata[[3]][[x]][, sim])
  }
  pick1 <- function(x, rdata, sim) {
    return(rdata[[4]][[x]][, sim])
  }
  alln_0 <- sapply(names(rdata[[3]]), pick0, rdata = rdata, sim = sim)
  alln_1 <- sapply(names(rdata[[4]]), pick1, rdata = rdata, sim = sim)
  rdata_sim <- list(n_0 = alln_0, n_1 = alln_1)
  nint <- length(allint)
  ## Get estimates of vbeta and beta using outcome data and each
  ## interim nint
  estvbeta <- sapply(1:nint, .getParams, rdatasim = rdata_sim, odatasim = odata[[sim]][3:4],
    s = s, nint = nint, vcovmat = vcovmat)
  colnames(estvbeta) <- allint
  return(t(estvbeta))
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------

#-------------------------------------------------------------------------------
# 7.  .getMeanvar
#-------------------------------------------------------------------------------
.getMeanvar <- function(mod, mean = NULL, cmodel = NULL, sd = NULL, rho = NULL) {

  ## Set-up
  tfu <- mod$rdata$tfu
  s <- mod$rdata$s

  ## Variance matrix Correlation matrix
  if (is.null(cmodel)) {
    ## if cmodel==NULL, then cannot reset rho
    corrmat <- mod$idata$cmodel$corrmat
  } else {
    ## Validate cmodel If rho=NULL, accept set value
    if (is.null(rho)) {
      rho <- mod$idata$cmodel$rho
    }
    ## Otherwise validate rho and reset to new value
    .valCorrmod(rmodel = "none", cmodel = cmodel, sd = mod$idata$sd,
      rho = rho, s = s)
    ## Get covariance matrix
    if (is.matrix(cmodel) == TRUE) {
      corrmat <- cmodel
    } else {
      if (cmodel == "uniform") {
        corrmat <- corrUnif(rho = rho, tfu = tfu)
      } else if (cmodel == "exponential") {
        corrmat <- corrExp(rho = rho, tfu = tfu)
      }
    }
  }
  ## SD matrix
  if (is.null(sd)) {
    ## If SD=NULL, accept set value
    sdmat <- diag(rep(mod$idata$sd, s))
  } else {
    ## Otherwise validate SD and reset to new value
    .valCorrmod(rmodel = "none", cmodel = mod$idata$cmodel$type, sd = sd,
      rho = mod$idata$cmodel$rho, s = s)
    if (length(sd) == 1) {
      sdmat <- diag(rep(sd, s))
    } else {
      sdmat <- diag(sd)
    }
  }
  vcovmat <- sdmat %*% corrmat %*% sdmat
  colnames(vcovmat) <- rownames(vcovmat) <- tfu$tfu
  colnames(corrmat) <- rownames(corrmat) <- tfu$tfu

  ## Mean matrix
  if (is.null(mean)) {
    meanmat_0 <- matrix(rep(0, s), nrow = 1, ncol = s)
    meanmat_1 <- matrix(c(rep(0, s - 1), mod$power$theta), nrow = 1,
      ncol = s)
  } else {
    ## validate nmean
    .valMeanmod(mean, s = s)
    meanmat_0 <- mean[1, ]
    meanmat_1 <- mean[2, ]
  }
  meanmat <- rbind(meanmat_0, meanmat_1)
  colnames(meanmat) <- tfu$tfu
  rownames(meanmat) <- names(mod$rdata$n)[1:2]

  ## Output
  return(list(vcovmat = vcovmat, corrmat = corrmat, meanmat = meanmat))
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------

#-------------------------------------------------------------------------------
# 8.  .getNinterim
#-------------------------------------------------------------------------------
.getNinterim <- function(x, rdat, tinterims, pick = NULL, nsim) {
  getInt <- function(y, x, rdat) {
    return(rdat[["n"]][[x]][rdat[["t"]] == y])
  }
  get_nint <- t(sapply(tinterims, getInt, x = x, rdat = rdat))
  if (nsim == 1) {
    get_nint <- t(get_nint)
  }
  if (!is.null(pick)) {
    get_nint <- as.matrix(get_nint[, !pick])
  }
  if (!all(is.na(get_nint))) {
    rownames(get_nint) <- tinterims
    colnames(get_nint) <- 1:dim(get_nint)[2]
  } else {
    get_nint <- NULL
  }
  return(get_nint)
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------

#-------------------------------------------------------------------------------
# 9.  .getParams
#-------------------------------------------------------------------------------
.getParams <- function(iint, rdatasim, odatasim, s, nint, vcovmat) {
  ## Get numbers
  n_0 <- rdatasim$n_0[iint, 2:(s + 1)]
  n_1 <- rdatasim$n_1[iint, 2:(s + 1)]
  ## Get totals for outcomes
  sums_0 <- odatasim[[1]][iint][[1]]
  sums_1 <- odatasim[[2]][iint][[1]]
  if (n_0[s] < 1 | n_1[s] < 1) {
    vbeta <- beta <- NA
  } else {
    n_0 <- c(n_0, 0)
    n_1 <- c(n_1, 0)
    beta_sum <- matrix(0, nrow = 2 * s, ncol = 1)
    V0 <- V1 <- matrix(0, nrow = s, ncol = s)
    for (r in 1:s) {
      for (c in 1:s) {
        mx <- max(r, c)
        for (mx in max(r, c):s) {
          iest_vcov <- solve(vcovmat[1:mx, 1:mx])
          V1[r, c] <- V1[r, c] + (n_1[mx] - n_1[mx + 1]) * iest_vcov[r,
          c]
          V0[r, c] <- V0[r, c] + (n_0[mx] - n_0[mx + 1]) * iest_vcov[r,
          c]
        }
      }
      xisigma <- rbind(solve(vcovmat[1:r, 1:r]), matrix(0, nrow = s -
        r, ncol = r))
      xisigma_0 <- rbind(xisigma, matrix(0, nrow = s, ncol = r))
      xisigma_1 <- rbind(xisigma, xisigma)
      beta_sum <- beta_sum + xisigma_0 %*% as.matrix(sums_0[1:r,
        r]) + xisigma_1 %*% as.matrix(sums_1[1:r, r])
    }
    iV0 <- solve(V0)
    iV1 <- solve(V1)
    iV01 <- iV0 + iV1
    ## Get values for beta and variance beta
    vbeta <- rbind(cbind(iV0, -iV0), cbind(-iV0, iV01))
    beta <- vbeta %*% beta_sum
  }
  ## Return values for treat effect
  return(c(vbeta = vbeta[2 * s, 2 * s], beta = beta[2 * s]))
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------

#-------------------------------------------------------------------------------
# 11.  gsearlyFit
#-------------------------------------------------------------------------------
gsearlyFit <- function (data,
          datanames = c("id", "atime", "catime", "intervention", "outcome"),
          cmodel=FALSE, vmodel=TRUE, method = "ML", full = FALSE){

  ## Construct model formula
  Call <- match.call()

  ## Validate inputs
  if (is.data.frame(data) == FALSE) {
    stop("data: must be a data.frame")
  }
  data <- data[, datanames]
  if (dim(data)[2] != 5) {
    stop("data: must have only the following five variables; id, atime, catime, 
                                intervention and outcome")
  }
  datanames <- names(data)
  s <- nlevels(data[, 2])
  add_term <- paste("(", datanames[2], "-1", ")", sep = "")
  int_term <- paste(paste("(", datanames[2], "-1", ")", sep = ""),
                    datanames[4], sep = ":")
  mformula <- as.formula(paste(datanames[5], "~", paste(add_term,
                                                        int_term, sep = "+")))
  cformula <- as.formula(paste("~ 1|", datanames[1], sep = ""))
  ccformula <- as.formula(paste("~ ", datanames[3], "|", datanames[1], sep = ""))
  vformula <- as.formula(paste("~ 1|", datanames[2], sep = ""))

  ## Fit gls
  ## Check method
  if (is.element(method, c("ML", "REML")) == FALSE) {
    stop("method: must be set to either ML or REML")
  }
  ## For selected correlation model; corSymm, corCompSym or corCAR1
  if(cmodel == "uniform"){
    if(vmodel == TRUE){
      gls_fit <- gls(mformula, data = data,
                      correlation = corCompSymm(form = cformula),
                      weights = varIdent(form = vformula), method = method)
    } else {
      gls_fit <- gls(mformula, data = data,
                           correlation = corCompSymm(form = cformula),
                           weights = NULL, method = method)
    }
  } else if (cmodel == "exponential"){
    if(vmodel == TRUE){
      gls_fit <- gls(mformula, data = data,
                      correlation = corCAR1(form = ccformula),
                      weights = varIdent(form = vformula), method = method)
    } else {
      gls_fit <- gls(mformula, data = data,
                           correlation = corCAR1(form = ccformula),
                           weights = NULL, method = method)
    }
  } else if (cmodel == FALSE){
    if(vmodel == TRUE){
      gls_fit <- gls(mformula, data = data,
                      correlation = corSymm(form = cformula),
                      weights = varIdent(form = vformula), method = method)
    } else {
      gls_fit <- gls(mformula, data = data,
                           correlation = corSymm(form = cformula),
                           weights = NULL, method = method)
    }
  }

  ## Extract vbeta and beta and correlation and covariance matrices
  if (full == FALSE) {
    N <- gls_fit$dims$N
    p <- gls_fit$dims$p
    s <- as.numeric(nlevels(data[,datanames[2]]))
    estcorr <- as.numeric(coef(gls_fit$model$corStruct, unconstrained = FALSE))
    estpsd <- as.numeric(coef(gls_fit$model$varStruct, unconstrained = FALSE))
    if(vmodel == TRUE){
      estsd <- c(1, as.numeric(estpsd)) * gls_fit$sigma
    } else {
      estsd <- gls_fit$sigma
    }
    evbeta <- as.numeric(vcov(gls_fit)[2 * s, 2 * s])
    ebeta <- as.numeric(coef(gls_fit)[2 * s])
    ez <- as.numeric(ebeta/sqrt(evbeta))
    outdat <- list(parameters = c(vbeta = evbeta, beta = ebeta,
        z = ez), model = list(N = N, p = p, estcorr = estcorr, estsd = estsd))
  }
  ## Output
  if (full == FALSE) {
    return(outdat)
  }
  else {
    gls_fit$call$data <- Call$data
    gls_fit$call$model <- get(as.character(gls_fit$call$model))
    gls_fit$call$correlation$form <- get(as.character(gls_fit$call$correlation$form))
    if(vmodel == TRUE){
      gls_fit$call$weights$form <- get(as.character(gls_fit$call$weights$form))
    }
    return(gls_fit)
  }
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------

#-------------------------------------------------------------------------------
# 15.  gsearlySimulate
#-------------------------------------------------------------------------------
gsearlySimulate <- function(mod, nsim = 1, minsamp = c(1,1), mean = NULL,
                     cmodel = NULL, sd = NULL, rho = NULL, digits = 3, 
                     odigits = 3, full = FALSE) {

  ## Simulation works only if recruit model is not set to 'none'
  if (mod$rdata$rmodel == "none") {
    stop("rmodel: can not simulate data if rmodel is set to none")
  }

  ## Need numbers recruited at interims, so set intonly==TRUE
  recruitn <- .simRdata(mod = mod, nsim = nsim, tinterims = mod$rdata$tinterims,
                        minsamp = minsamp, intonly = TRUE)

  ## Mean and covariance matrices, and rho and sd for simulating
  ## outcomes
  meanvar <- .getMeanvar(mod = mod, mean = mean, cmodel = cmodel, sd = sd,
                         rho = rho)

  ## Outcome data
  odata <- lapply(1:nsim, .simOdata, mod = mod, recruitn = recruitn,
                  meanmat = meanvar$meanmat, vcovmat = meanvar$vcovmat,
                  odigits = odigits)
  if (full == TRUE) {
    data_0 <- lapply(odata, "[[", 1)
    data_1 <- lapply(odata, "[[", 2)
  }

  ## Estimate varB and B
  allint <- recruitn$tinterims
  nint <- length(recruitn$tinterims)
  param_est <- lapply(1:nsim, .estParams, rdata = recruitn, odata = odata,
                      s = mod$rdata$s, vcovmat = meanvar$vcovmat, allint = allint)
  fnVbeta <- function(x, param_est) {
    return(unlist(lapply(param_est, "[[", x)))
  }
  vbeta <- sapply(1:nint, fnVbeta, param_est = param_est)
  beta <- sapply((nint + 1):(nint + nint), fnVbeta, param_est = param_est)
  beta <- matrix(as.numeric(beta), nrow=nsim, ncol=mod$rdata$s, byrow=FALSE)
  vbeta <- matrix(as.numeric(vbeta), nrow=nsim, ncol=mod$rdata$s, byrow=FALSE)
  colnames(vbeta) <- colnames(beta) <- allint
  rownames(vbeta) <- rownames(beta) <- 1:nsim

  ## Output
  simmodel <- list(nsim = recruitn$nsim, s = mod$rdata$s, tlooks = allint,
                   tfu = mod$rdata$tfu, mean = meanvar$meanmat, covariance = meanvar$vcovmat,
                   correlation = meanvar$corrmat)
  if (full == TRUE) {
    nout <- list(recruitn[[3]], recruitn[[4]])
    dataout <- list(data_0, data_1)
    names(nout) <- names(dataout) <- names(mod$rdata$n)[1:2]
    outdat <- list(model = simmodel, full = full, 
                   parameters = list(vbeta = round(vbeta, digits = digits),
                             beta = round(beta, digits = digits)), n = nout, data = dataout)
    return(outdat)
  } else {
    return(list(model = simmodel, full = full, 
                   parameters = list(vbeta = round(vbeta, digits = digits),
                             beta = round(beta, digits = digits))))
  }
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------

#-------------------------------------------------------------------------------
# 17.  modelParameters
#-------------------------------------------------------------------------------
modelParameters <- function (data,
      datanames = c("id", "atime", "intervention", "outcome"), vcovmat) {
  
  ## Validate inputs
  if (is.data.frame(data) == FALSE) {
    stop("data: must be a data.frame")
  }
  data <- data[, datanames]
  if (dim(data)[2] != 4) {
    stop("data: must have only the following four variables; id, atime, intervention and outcome")
  }
  datanames <- names(data)
  id <- data[, datanames[1]]
  s <- nlevels(data[, 2])
  if (!isSymmetric(vcovmat) | any(eigen(vcovmat)$values < 0)) {
    stop("vcovmat: must be a valid covariance matrix")
  }
  if (dim(vcovmat)[1] != s | dim(vcovmat)[2] != s) {
    stop("vcovmat: set to s x s matrix")
  }
  add_term <- paste("(", datanames[2], "-1", ")", sep = "")
  int_term <- paste(paste("(", datanames[2], "-1", ")", sep = ""),
                    datanames[3], sep = ":")
  mformula <- as.formula(paste(datanames[4], "~", paste(add_term,
                                                        int_term, sep = "+")))
  dmatrix <- model.matrix(mformula, data = data)
  mod_mat <- data.frame(id = as.numeric(id), dmatrix)
  split_dmatrix <- split(as.data.frame(mod_mat[, 2:((2 * s) +
                                                      1)]), mod_mat$id)
  n <- length(split_dmatrix)
  split_data <- split(as.data.frame(data), id)
  beta_sum <- matrix(0, nrow = 2 * s, ncol = 1)
  var_sum <- matrix(0, nrow = 2 * s, ncol = 2 * s)
  for (i in 1:n) {
    nout <- dim(split_dmatrix[[i]])[1]
    var_sum <- var_sum + t(as.matrix(split_dmatrix[[i]])) %*%
      solve(vcovmat[1:nout, 1:nout]) %*% as.matrix(split_dmatrix[[i]])
    beta_sum <- beta_sum + t(as.matrix(split_dmatrix[[i]])) %*%
      solve(vcovmat[1:nout, 1:nout]) %*% as.matrix(split_data[[i]][,
                                                                   datanames[4]])
  }
  estvar <- solve(var_sum)
  estbeta <- estvar %*% beta_sum
  evbeta <- as.numeric(estvar[2 * s, 2 * s])
  ebeta <- as.numeric(estbeta[2 * s])
  ez <- as.numeric(ebeta/sqrt(evbeta))
  outdat <- c(vbeta = evbeta, beta = ebeta, z = ez)
  return(outdat)
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------

#-------------------------------------------------------------------------------
# 30.  simdataExtract
#-------------------------------------------------------------------------------
simdataExtract <- function(simmod, simn, tlooks, datanames = c("id",
  "atime", "catime", "intervention", "outcome"), full = FALSE) {

  ## Validate inputs
  if(simmod$full == FALSE){
    stop("full: set full to TRUE in the simulation model")   
  }
  if (is.element(simn, 1:simmod$model$nsim) == FALSE) {
    stop("simn: invalid simulation number")
  }
  if (is.element(tlooks, simmod$model$tlooks) == FALSE) {
    stop("tlooks: invalid interim time")
  }

  ## Pick simulation number
  get_int <- which(simmod$model$tlooks == tlooks)
  n0 <- max(simmod$n[[1]][[1]][, simn], na.rm = TRUE)
  n1 <- max(simmod$n[[2]][[1]][, simn], na.rm = TRUE)
  s <- length(simmod$model$tfu$tfu)

  ## Outcomes Group 0
  outcome_0 <- simmod$data[[1]][[simn]][get_int]
  atime_0 <- factor(rep(simmod$model$tfu$tfu, each = n0),ordered=TRUE)
  catime_0 <- rep(simmod$model$tfu$stfu, each = n0)
  id_0 <- rep(1:n0, times = s)
  intervention_0 <- rep(names(simmod[["data"]])[1], s * n0)
  simdat_0 <- data.frame(id = id_0, atime = atime_0, catime = catime_0,
    intervention = intervention_0, outcome = as.numeric(outcome_0[[1]]))
  simdat_0 <- subset(simdat_0, complete.cases(simdat_0))
  ## Group 1
  outcome_1 <- simmod$data[[2]][[simn]][get_int]
  atime_1 <- factor(rep(simmod$model$tfu$tfu, each = n1),ordered=TRUE)
  catime_1 <- rep(simmod$model$tfu$stfu, each = n1)
  id_1 <- rep(max(simdat_0$id) + 1:n1, times = s)
  intervention_1 <- rep(names(simmod[["data"]])[2], s * n1)
  simdat_1 <- data.frame(id = id_1, atime = atime_1, catime = catime_1,
    intervention = intervention_1, outcome = as.numeric(outcome_1[[1]]))
  simdat_1 <- subset(simdat_1, complete.cases(simdat_1))

  ## Bind data together and re-order
  simdat <- rbind(simdat_0, simdat_1)
  simdat$intervention <- factor(simdat$intervention,
                          levels=rownames(simmod$model$mean))
  simdat <- simdat[order(simdat$id), ]

  ## Get meta data
  if (full == TRUE) {
    vbeta <- as.numeric(simmod$parameters$vbeta[simn, get_int])
    beta <- as.numeric(simmod$parameters$beta[simn, get_int])
    model <- simmod$model[-1]
    model$tlooks <- simmod$model$tlooks[get_int]
    model$parameters <- c(vbeta = vbeta, beta = beta)
    model$n <- t(colSums(table(simdat$id, simdat$atime, simdat$intervention)))
  }

  ## Output
  rownames(simdat) <- NULL
  colnames(simdat) <- datanames
  if (full == TRUE) {
    outdat <- list(model = model, data = simdat)
  } else {
    outdat <- simdat
  }
  return(outdat)
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------

#-------------------------------------------------------------------------------
# 31.  .simOdata
#-------------------------------------------------------------------------------
.simOdata <- function(x, mod, recruitn, meanmat, vcovmat, odigits) {

  ## Set-up
  s <- mod$rdata$s
  tfu <- mod$rdata$tfu
  allint <- recruitn$tinterims
  nint <- length(allint)

  ## Numbers recruited
  n_0 <- recruitn[3][[1]][[1]][nint, x]
  n_1 <- recruitn[4][[1]][[1]][nint, x]

  ## Outcomes
  outcomes_0 <- round(rmvnorm(as.integer(n_0), mean = meanmat[1, ],
    sigma = vcovmat), odigits)
  outcomes_1 <- round(rmvnorm(as.integer(n_1), mean = meanmat[2, ],
    sigma = vcovmat), odigits)


  ## Functions for data summaries
  getSums <- function(q, y, s, dat) {
    sumvec <- c(colSums(matrix(dat[1:y[q], 1:q], ncol = q, byrow = FALSE)),
      rep(0, s - q))
    return(sumvec)
  }
  diffSum <- function(q, mat, s) {
    if (q == s) {
      mat[, q] <- mat[, q]
    } else {
      mat[, q] <- mat[, q] - mat[, q + 1]
    }
    names(mat[, q]) <- 1:s
    return(mat[, q])
  }
  getN <- function(q, dat, int, ss) {
    return(dat[[q + 1]][int, ss])
  }
  makeNA <- function(q, nj, dat) {
    new_dat <- dat
    if (nj[q] < dim(dat)[1]) {
      new_dat[(nj[q] + 1):dim(dat)[1], q] <- NA
    }
    return(new_dat[, q])
  }
  setZero <- function(mat, s) {
    new_mat <- matrix(0, ncol = s, nrow = s)
    new_mat[s, s] <- mat[s, s]
    return(new_mat)
  }

  ## Get sums at all interims and create missingness patterns
  tot_sum_0 <- tot_sum_1 <- vector("list", nint)
  new_data_0 <- new_data_1 <- vector("list", nint)
  names(tot_sum_0) <- names(tot_sum_1) <- names(new_data_0) <- names(new_data_1) <- allint
  for (j in nint:1) {
    nj_0 <- sapply(1:s, getN, dat = recruitn[3][[1]], int = j, ss = x)
    nj_1 <- sapply(1:s, getN, dat = recruitn[4][[1]], int = j, ss = x)
    new_data_0[[j]] <- sapply(1:s, makeNA, nj = nj_0, dat = outcomes_0)
    new_data_1[[j]] <- sapply(1:s, makeNA, nj = nj_1, dat = outcomes_1)
    tot_sum_0[[j]] <- sapply(1:s, getSums, y = nj_0, s = s, dat = new_data_0[[j]])
    tot_sum_1[[j]] <- sapply(1:s, getSums, y = nj_1, s = s, dat = new_data_1[[j]])
    if (j == nint) {
      tot_sum_0[[j]] <- setZero(mat = tot_sum_0[[j]], s = s)
      tot_sum_1[[j]] <- setZero(mat = tot_sum_1[[j]], s = s)
    }
    ## Create matrices of sums
    tot_sum_0[[j]][lower.tri(tot_sum_0[[j]])] <- NA
    tot_sum_1[[j]][lower.tri(tot_sum_1[[j]])] <- NA
    tot_sum_0[[j]] <- sapply(1:s, diffSum, mat = tot_sum_0[[j]], s = s)
    tot_sum_1[[j]] <- sapply(1:s, diffSum, mat = tot_sum_1[[j]], s = s)
    rownames(new_data_0[[j]]) <- 1:n_0
    rownames(new_data_1[[j]]) <- 1:n_1
    colnames(new_data_0[[j]]) <- colnames(new_data_1[[j]]) <- as.character(mod$rdata$tfu$tfu)
    rownames(tot_sum_0[[j]]) <- rownames(tot_sum_1[[j]]) <- as.character(mod$rdata$tfu$tfu)
    colnames(tot_sum_0[[j]]) <- colnames(tot_sum_1[[j]]) <- as.character(mod$rdata$tfu$tfu)
  }

  ## Output
  odata <- list(new_data_0, new_data_1)
  ototals <- list(tot_sum_0, tot_sum_1)
  names(odata) <- names(ototals) <- names(mod$rdata$n)[1:2]
  outdat <- c(data = odata, totals = ototals)
  return(outdat)
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------

#-------------------------------------------------------------------------------
# 32.  .simRdata
#-------------------------------------------------------------------------------
.simRdata <- function(mod, tinterims = NULL, nsim = 1, minsamp = c(1, 1),
  intonly = TRUE) {

  ## Set-up
  getrmodel <- .selectRmodel(rmodel = mod$rdata$rmodel)
  n <- as.numeric(mod$rdata$n[3])
  n_0 <- as.numeric(mod$rdata$n[1])
  n_1 <- as.numeric(mod$rdata$n[2])
  trecruit <- mod$rdata$trecruit
  tfu <- mod$rdata$tfu$tfu
  alltfu <- c(0, tfu)
  s <- mod$rdata$s
  ss <- s + 1
  m <- mod$rdata$m

  ## Check tinterims
  if (is.null(tinterims)) {
    tinterims <- mod$rdata$tinterims
  } else {
    .valTinterims(tinterims)
    i_tinterim <- any(tinterims <= tfu[s]) || any(tinterims >= trecruit)
    if (i_tinterim == TRUE) {
      stop("tinterims: infeasible tinterims, need tfu[s] < tinterims < trecruit")
    }
  }
  allint <- c(tinterims, trecruit + tfu[s])

  ## Get fixed or random numbers, based on poisson
  rate_int_0 <- sapply(0:(trecruit + tfu[s]), getrmodel$fr, n = n_0,
    trecruit = trecruit, m = m)
  rate_int_1 <- sapply(0:(trecruit + tfu[s]), getrmodel$fr, n = n_1,
    trecruit = trecruit, m = m)
  sim_recruit_0 <- sapply(rate_int_0, rpois, n = nsim)
  sim_recruit_1 <- sapply(rate_int_1, rpois, n = nsim)

  ## Cumulative recruitment numbers
  if (nsim == 1) {
    csim_recruit_0 <- matrix(cumsum(sim_recruit_0), ncol = 1)
    csim_recruit_1 <- matrix(cumsum(sim_recruit_1), ncol = 1)
  } else {
    mcumsum <- function(x, dat) {
      return(cumsum(dat[x, ]))
    }
    csim_recruit_0 <- sapply(1:nsim, mcumsum, dat = sim_recruit_0)
    csim_recruit_1 <- sapply(1:nsim, mcumsum, dat = sim_recruit_1)
  }
  colnames(csim_recruit_0) <- colnames(csim_recruit_1) <- 1:nsim
  rownames(csim_recruit_0) <- rownames(csim_recruit_1) <- 0:(trecruit +
    tfu[s])

  ## Get follow-up numbers from recruited numbers
  ffu <- function(x, tfu, s, nrecruit, trecruit, nsim) {
    data_mat <- rbind(matrix(0, nrow = tfu[x], ncol = nsim), as.matrix(nrecruit[1:(trecruit +
      1 + tfu[s] - tfu[x]), ]))
    colnames(data_mat) <- 1:nsim
    rownames(data_mat) <- 0:(trecruit + tfu[s])
    return(data_mat)
  }
  rdat_0 <- lapply(1:ss, ffu, tfu = alltfu, s = ss, nrecruit = csim_recruit_0,
    trecruit = trecruit, nsim = nsim)
  rdat_0 <- list(t = 0:(trecruit + tfu[s]), n = rdat_0)
  rdat_1 <- lapply(1:ss, ffu, tfu = alltfu, s = ss, nrecruit = csim_recruit_1,
    trecruit = trecruit, nsim = nsim)
  rdat_1 <- list(t = 0:(trecruit + tfu[s]), n = rdat_1)
  fu_names <- as.character(alltfu)
  names(rdat_0[["n"]]) <- names(rdat_1[["n"]]) <- fu_names

  ## Get interim numbers Exclude simulations where numbers at first
  ## interim are less than minsamp
  testmin_0 <- .getNinterim(fu_names[ss], rdat = rdat_0, tinterims = allint,
    nsim = nsim)
  testmin_1 <- .getNinterim(fu_names[ss], rdat = rdat_1, tinterims = allint,
    nsim = nsim)
  ipick <- as.logical(as.numeric(testmin_0[1, ] < minsamp[1]) + as.numeric(testmin_1[1,
    ] < minsamp[2]))
  getinterims_0 <- lapply(fu_names, .getNinterim, rdat = rdat_0, tinterims = allint,
    pick = ipick, nsim = nsim)
  getinterims_1 <- lapply(fu_names, .getNinterim, rdat = rdat_1, tinterims = allint,
    pick = ipick, nsim = nsim)
  names(getinterims_0) <- names(getinterims_1) <- fu_names
  nnsim <- nsim - sum(ipick)

  ## Output
  if (intonly == FALSE) {
    simdat <- list(rmodel = mod$rdata$rmodel, trecruit = trecruit,
      s = s, tfu = mod$rdata$tfu, n = n, vphi = mod$rdata$vphi, m = m,
      nsim = nnsim, data = list(rdat_0, rdat_1), tinterims = tinterims,
      interims = list(getinterims_0, getinterims_1))
    names(simdat[["data"]]) <- names(simdat[["interims"]]) <- names(mod$rdata$n)[1:2]
  } else {
    simdat <- list(nsim = nnsim, tinterims = allint, getinterims_0,
      getinterims_1)
    names(simdat)[3:4] <- names(mod$rdata$n)[1:2]
  }
  return(simdat)
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------

#-------------------------------------------------------------------------------
# 49.  dataOrder
#-------------------------------------------------------------------------------
dataOrder <- function(data,
          datanames = c("id", "atime", "catime", "intervention", "outcome")){
  ## Sort data into appropriate form
  id_list <- split(data, data[, datanames[1]])
  nlist <- length(id_list)
  ordersubj <- function(x, id_list, datanames) {
    id_list[[x]][order(id_list[[x]][, datanames[2]]), ]
  }
  lodata <- lapply(1:nlist, ordersubj, id_list = id_list, datanames = datanames)
  odata <- Reduce(rbind, lodata)
  return(odata)
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------

#-------------------------------------------------------------------------------
# end
#-------------------------------------------------------------------------------