R/TwoSample.Wald.and.Boundary.R

In gsMeanFreq: Group Sequential Clinical Trial Designs for Composite Endpoints

#' @title Function to calculate the Wald statistics and group sequential boundaries for two-sample monitoring.
#' @description
#' Computes stage-wise Wald statistics from estimated two-sample statistics and their variances, determines information
#' fractions, allocates Type I error using a user specified spending function, and constructs corresponding two-sided
#' rejection boundaries. The function supports both a consistent correlation approach (supplied in \code{Q.cov.est}) and
#' a canonical correlation approach based on information fractions.
#'
#'
#' @param Q.cov.est A list containing stage-wise estimates and covariances, either returned by \code{TwoSample.Estimator.LR.sequential()}
#' or \code{TwoSample.Estimator.GT.sequential()}.
#' @param spend A function specifying the cumulative Type I error spending function \eqn{\alpha(t)} evaluated at information fraction
#' \code{t}. For example, \code{OBF}.
#' @param calendars Numeric vector of analysis calendar times (in years), defining the planned monitoring schedule and number of analysis.
#' @param alpha Overall two-sided Type I error.
#' @param planned.n Planned total sample size at the final analysis.
#' @param Iunit Information per subject (or per unit of sample size) used to scale the information fractions.
#'
#' @returns A list with components:
#' \itemize{
#' \item \code{Qs}: Stage-wise test statistics.
#' \item \code{vars}: Stage-wise variance estimates for \code{Qs}.
#' \item \code{raw.information}: Information fractions prior to any adjustments for early completion or skipped analysis.
#' \item \code{Wald}: Stage-wise Wald statistics \code{Qs/sqrt(vars)}.
#' \item \code{consistent.bdry}: Two-sided rejection boundaries computed using the consistent correlation matrix.
#' \item \code{canonical.bdry}: Two-sided rejection boundaries computed using the canonical correlation matrix.
#' \item \code{consistent.reject}: Indicator vector for boundary crossing under the consistent approach (only the the first crossing is retained).
#' \item \code{canonical.reject}: Indicator vector for boundary crossing under the canonical approach (only the the first crossing is retained).
#' \item \code{nu}: Final information fractions used for boundary construction.
#' \item \code{pi}: Incremental Type I error allocated to each analysis.
#' \item \code{total.ns}: Total accrued sample size at each analysis.
#' }
#' @export
#' @importFrom stats qnorm
#'
TwoSample.Wald.and.Boundary <- function(Q.cov.est, spend, calendars, alpha, planned.n, Iunit){

  Qs <- Q.cov.est$Qs
  vars <- Q.cov.est$vars

  consistent.corr <- Q.cov.est$corr.matrix # the consistent method's correlation matrix is already calculated

  # nss <- Q.cov.est$nss
  total.ns <- Q.cov.est$total.ns
  nu <- c()
  pi <- c()
  consistent.bdry <- c()
  canonical.bdry <- c()


  ##### Calculate the Wald statistics #####
  Wald <- Qs/sqrt(vars)
  # transformed.Wald <- log(Qs) - log(sqrt(vars))

  # Step 1: Calculate the information fraction for this simulated data
  for (j in 1:(length(calendars))){
    nu[j] = vars[j]/(planned.n*Iunit)
    nu[length(calendars)] = 1
  }

  # Qinghua 04/22/2025 Update: calculate the canonical correlation matrix
  raw.information <- nu
  canonical.corr <- diag(3)
  for (p in 1:(3-1)){
    for (q in (p+1):3) {
      canonical.corr[p,q] <- sqrt(raw.information[p]/raw.information[q])
      canonical.corr[q,p] <- canonical.corr[p,q]
    }
  }

  # Step 2: Calculate the allocation of type I error and rejection boundaries
  if (nu[1] >= 1){
    # Reached 100% information at stage 1, in that case, treat the first analysis as the final one.
    # Carried the 100% information fraction forward and set the stages 2 and 3 boundaries to 99
    nu = c(1, 1, 1)
    pi = c(alpha, 0, 0)
    consistent.bdry = c(qnorm(p = alpha/2, lower.tail = FALSE), 99, 99)
    canonical.bdry = c(qnorm(p = alpha/2, lower.tail = FALSE), 99, 99)
    Wald = c(Wald[1], 0, 0)
    # transformed.Wald = c(transformed.Wald[1], 0, 0)

  } else if (nu[2] >= 1){
    # Reached 100% information at stage 2, in that case, treat the 2nd interim analysis as the final one.
    # Carried the 100% information fraction forward and set the stage 3 boundary to 99
    nu = c(nu[1], 1, 1)
    pi = c(spend(nu[1]), alpha - spend(nu[1]), 0)

    consistent.bdry[1] = qnorm(p = pi[1]/2, lower.tail = FALSE)
    consistent.bdry[2] = find.c(previous.c = consistent.bdry[1], pi = pi[2], corr = consistent.corr[c(1:2),c(1:2)])
    consistent.bdry[3] = 99

    canonical.bdry[1] = qnorm(p = pi[1]/2, lower.tail = FALSE)
    canonical.bdry[2] = find.c(previous.c = canonical.bdry[1], pi = pi[2], corr = canonical.corr[c(1:2),c(1:2)])
    canonical.bdry[3] = 99

    Wald = c(Wald[1], Wald[2], 0)
    # transformed.Wald = c(transformed.Wald[1], transformed.Wald[1], 0)

  } else if (nu[2]< nu[1]){
    # Observed a smaller information fraction at a later stage
    # in that case, skip stage 2 and move on to the final stage
    # set the stage 2 boundary to 88
    nu = c(nu[1], nu[1], 1)
    pi = c(spend(nu[1]), 0, alpha - spend(nu[1]))

    consistent.bdry[1] = qnorm(p = pi[1]/2, lower.tail = FALSE)
    consistent.bdry[2] = 88
    consistent.bdry[3] = find.c(previous.c = consistent.bdry[1], pi = pi[3], corr = consistent.corr[c(1,3),c(1,3)])

    canonical.bdry[1] = qnorm(p = pi[1]/2, lower.tail = FALSE)
    canonical.bdry[2] = 88
    canonical.bdry[3] = find.c(previous.c = canonical.bdry[1], pi = pi[3], corr = canonical.corr[c(1,3),c(1,3)])

    Wald = c(Wald[1], 0, Wald[3])
    # transformed.Wald = c(transformed.Wald[1], 0, transformed.Wald[3])
  } else {
    # Trial proceed as planned
    # Calculate the allocation of Type I errors
    pi[1] <- spend(nu[1])
    for (j in 2:(length(calendars) - 1)){
      pi[j] <- spend(nu[j]) - spend(nu[j-1])
    }
    pi[length(calendars)] <- alpha - spend(nu[length(calendars) - 1])

    #  Calculate the rejection boundary
    consistent.bdry[1] = qnorm(p = pi[1]/2, lower.tail = FALSE)
    for (j in 2:length(calendars)){
      consistent.bdry[j] <- find.c(previous.c = consistent.bdry[1:j-1], pi = pi[j], corr = consistent.corr[1:j,1:j])
    }

    canonical.bdry[1] = qnorm(p = pi[1]/2, lower.tail = FALSE)
    for (j in 2:3){
      canonical.bdry[j] <- find.c(previous.c = canonical.bdry[1:j-1], pi = pi[j], corr = canonical.corr[1:j,1:j])
    }
  }

  # Step 4: Compare the Wald statistics and the rejection boundary
  consistent.reject <- 1*(Wald >= consistent.bdry|Wald <= -consistent.bdry)
  # transformed.reject = 1*(transformed.Wald >= bdry |transformed.Wald <= -bdry)
  canonical.reject <- 1*(Wald >= canonical.bdry|Wald <= -canonical.bdry)


  # Step 5: find the first reject and set the later values to zero
  if (any(consistent.reject == 1)) {
    first <- which(consistent.reject == 1)[1]
    if (first < length(consistent.reject)) {
      consistent.reject[(first + 1):length(consistent.reject)] <- 0
    }
  }

  if (any(canonical.reject == 1)){
    first <- which(canonical.reject == 1)[1]
    if (first < length(canonical.reject)){
      canonical.reject[(first + 1):length(canonical.reject)] <- 0
    }
  }

  return(list(Qs = Qs,
              vars = vars,
              raw.information = raw.information,
              Wald = Wald,
              # transformed.Wald = transformed.Wald,
              consistent.bdry = consistent.bdry,
              canonical.bdry = canonical.bdry,
              consistent.reject = consistent.reject,
              canonical.reject = canonical.reject,
              # transformed.reject = transformed.reject,
              nu = nu,
              pi = pi,
              total.ns = total.ns))
}