R/Group_specific_Delta_AUC_estimation.R

In marie-alexandre/AUCcomparison: Statistical test - Between group comparison of AUC

#' @title Difference of AUC of Two Group-Specific Polynomial Marginal Dynamics
#' @description \loadmathjax This function estimates the difference of area under the curve of marginal dynamics from two groups when marginal dynamics are modeled by group-structured polynomials or B-spline curves.

#' @param MEM_Pol_group A list with similar structure than the output provided by the function \link[AUCcomparison]{MEM_Polynomial_Group_structure}. 
#' 
#' A list containing: 
#' \itemize{
#' \item \code{Model_estimation}: a list containing at least 2 elements: \enumerate{
#'     \item the vector of the marginal (fixed) parameters estimates (at least for the groups whose AUC is to estimate), labeled _'beta'_.
#'     \item the variance-covariance matrix of these parameters, labeled _'varFix'_ (see \link[AUCcomparison]{MEM_Polynomial_Group_structure} for details about the parameter order).
#'     }
#' \item \code{Model_features}: a list of at least 2 elements: \enumerate{
#'     \item \code{Groups}: a vector indicating the names of the groups whose fixed parameters are given.
#'     \item \code{Marginal.dyn.feature}: a list summarizing the features of the marginal dynamics defined in the model: 
#'     \itemize{
#'         \item \code{dynamic.type}: a character scalar indicating the chosen type of marginal dynamics. Options are 'polynomial' or 'spline'.
#'         \item \code{intercept}: a logical vector summarizing choices about global and group-specific intercepts (Number of groups + 1) elements whose elements are named as ('global.intercept','group.intercept1', ..., 'group.interceptG') if G Groups are defined in \code{MEM_Pol_group}. For each element of the vector, if TRUE, the considered intercept is considered as included in the model (see \emph{Examples}).
#'         }
#'      If \code{dynamic.type} is defined as 'polynomial':\itemize{
#'          \item \code{polynomial.degree}: an integer vector indicating the degree of polynomial functions, one value for each group.
#'          }
#'      If \code{dynamic.type} is defined as 'spline':\itemize{
#'          \item \code{spline.degree}: an integer vector indicating the degree of B-spline curves, one for each group. 
#'          \item \code{knots}: a list of group-specific internal knots used to build B-spline basis (one numerical vector for each group) (see \link[splines]{bs} for more details).
#'          \item \code{df}: a numerical vector of group-specific degrees of freedom used to build B-spline basis, (one for each group).
#'          \item \code{boundary.knots}: a list of group-specific boundary knots used to build B-spline  basis (one vector for each group) (see \link[splines]{bs} for more details).
#'          }
#'     }
#' }
#' @param Group1 a character scalar indicating the name of the first group whose marginal dynamics must be considered. This group name must belong to the set of groups involved in the MEM (see \code{Groups} vector in \code{MEM_Pol_group}). 
#' @param Group2 a character scalar indicating the name of the second group whose marginal dynamics must be considered. This group name must belong to the set of groups involved in the MEM (see \code{Groups} vector in \code{MEM_Pol_group}).
#' @param time.G1 a numerical vector of time points (x-axis coordinates) to use for the Group1 AUC calculation.
#' @param time.G2 a numerical vector of time points (x-axis coordinates) to use for the Group2 AUC calculation.
#' @param common.interval a logical scalar. If FALSE, the difference of AUC is calculated as the difference of AUCs where the AUC of each group is calculated on its specific interval of time. If TRUE (default), the difference of AUC is estimated on a common interval of time defined as the intersect of the two group-specific interval (see @details for more details).
#' @param method a character scalar indicating the interpolation method to use to estimate the AUC. Options are 'trapezoid' (default), 'lagrange' and 'spline'. In this version, the 'spline' interpolation is implemented with "not-a-knot" spline boundary conditions.
#' @param Averaged a logical scalar. If TRUE, the function return the difference of normalized AUC (nAUC) where nAUC is computated as the AUC divided by the range of time of calculation. If FALSE (default), the classic AUC is calculated.
#' 
#' @details The difference of area under the curve between the two groups of interest is calculated as an approximation of the difference of the integrals of the expected value of the estimated outcome Y specific to the two groups g1 and g2.
#' Assuming a time interval \mjteqn{\[0,T_{g1}\]}{\[0,T_{g1}\]}{\[0,T_{g1}\]} for the group 1 and \mjteqn{\[0,T_{g2}\]}{\[0,T_{g2}\]}{\[0,T_{g2}\]} for the group 2, the difference of AUC is then calculated as
#' \mjtdeqn{\Delta AUC = \int_0^{T_{g_2}} E(\hat{Y_{g_2}})(t) dt - \int_0^{T_{g_1}} E(\hat{Y_{g_1}})(t) dt}{\Delta AUC = \int_0^{T_{g_2}} E(\hat{Y_{g_2}})(t) dt - \int_0^{T_{g_1}} E(\hat{Y_{g_1}})(t) dt}{\Delta AUC = \int_0^{T_{g_2}} E(\hat{Y_{g_2}})(t) dt - \int_0^{T_{g_1}} E(\hat{Y_{g_1}})(t) dt}
#' Similarly, the difference of the normalized AUC (\mjteqn{\Delta nAUC}{\Delta nAUC}{\Delta nAUC}) for these two same groups is then defined as 
#' \mjtdeqn{\Delta nAUC = \frac{1}{T_{g_2}}\int_0^{T_{g_2}} E(\hat{Y_{g_2}})(t) dt - \frac{1}{T_{g_1}}\int_0^{T_{g_1}} E(\hat{Y_{g_1}})(t) dt}{\Delta nAUC = \frac{1}{T_{g_2}}\int_0^{T_{g_2}} E(\hat{Y_{g_2}})(t) dt - \frac{1}{T_{g_1}}\int_0^{T_{g_1}} E(\hat{Y_{g_1}})(t) dt}{\Delta nAUC = \frac{1}{T_{g_2}}\int_0^{T_{g_2}} E(\hat{Y_{g_2}})(t) dt - \frac{1}{T_{g_1}}\int_0^{T_{g_1}} E(\hat{Y_{g_1}})(t) dt}
#' When we want to estimate the difference of AUC (or nAUC) on the same time interval for the two groups (\code{common.interval} = TRUE), we define \mjteqn{T = \min(T_{g_1},T_{g_2}))}{T = \min(T_{g_1},T_{g_2}))}{T = \min(T_{g_1},T_{g_2}))} and we calculate the difference of AUC as 
#' \mjtdeqn{\Delta AUC = \int_0^{T} E(\hat{Y_{g_2}})(t) dt - \int_0^{T} E(\hat{Y_{g_1}})(t) dt}{\Delta AUC = \int_0^{T} E(\hat{Y_{g_2}})(t) dt - \int_0^{T} E(\hat{Y_{g_1}})(t) dt}{\Delta AUC = \int_0^{T} E(\hat{Y_{g_2}})(t) dt - \int_0^{T} E(\hat{Y_{g_1}})(t) dt}
#' \mjtdeqn{\Delta nAUC = \frac{1}{T}\int_0^{T} E(\hat{Y_{g_2}})(t) dt - \frac{1}{T}\int_0^{T} E(\hat{Y_{g_1}})(t) dt}{\Delta nAUC = \frac{1}{T}\int_0^{T} E(\hat{Y_{g_2}})(t) dt - \frac{1}{T}\int_0^{T} E(\hat{Y_{g_1}})(t) dt}{\Delta nAUC = \frac{1}{T}\int_0^{T} E(\hat{Y_{g_2}})(t) dt - \frac{1}{T}\int_0^{T} E(\hat{Y_{g_1}})(t) dt}
#' 
#' @return A numerical scalar defined as \mjteqn{\Delta AUC = AUC_2 - AUC_1}{\Delta AUC = AUC_2 - AUC_1}{\Delta AUC = AUC_2 - AUC_1} (or \mjteqn{\Delta nAUC = nAUC_2 - nAUC_1}{\Delta nAUC = nAUC_2 - nAUC_1}{\Delta nAUC = nAUC_2 - nAUC_1})  with \mjteqn{AUC_1}{AUC_1}{AUC_1} (or \mjteqn{nAUC_1}{nAUC_1}{nAUC_1}) and  \mjteqn{AUC_2}{AUC_2}{AUC_2} (or \mjteqn{nAUC_2}{nAUC_2}{nAUC_2}) being respectively estimated as the AUC (or nAUC) for the Group1 and for the Group2.
#' 
#' @examples 
#' \donttest{# Download of data
#' data("HIV_Simu_Dataset_Delta01_cens")
#' data <- HIV_Simu_Dataset_Delta01_cens
#' 
#' # Change factors in character vectors
#' data$id <- as.character(data$id) ; data$Group <- as.character(data$Group)
#' 
#' # Example 1: We consider the variable \code{MEM_Pol_Group} as the output of our function 
#' # \link[AUCcomparison]{MEM_Polynomial_Group_structure}
#' MEM_estimation_1 <- MEM_Polynomial_Group_structure(y=data$VL,x=data$time,Group=data$Group,
#'                                                    Id=data$id,Cens=data$cens)
#' time_group1 <- unique(data$time[which(data$Group=="Group1")])
#' time_group2 <- unique(data$time[which(data$Group=="Group2")])
#'
#' Delta_AUC_1 <- Group_specific_Delta_AUC_estimation(MEM_Pol_group=MEM_estimation_1,
#'                                                    Group1="Group1",Group2="Group2",
#'                                                    time.G1=time_group1,time.G2=time_group2)
#'
#' # Example 2: We consider results of MEM estimation from another source. 
#' # We have to give build the variable 'MEM_Pol_group' with the good structure
#' # We build the variable 'MEM_estimation_2' with the results of MEM estimation obtained for 2 groups 
#' MEM_estimation_2 <- list(Model_estimation=c(1.077,0.858,-0.061,0.0013,0.887,-0.066,0.0014),
#'                         Model_features=list(Groups=c("Group1","Group2"),
#'                                             Marginal.dyn.feature=list(dynamic.type="polynomial",
#'                                                    intercept=c(global.intercept=TRUE,
#'                                                    group.intercept1=FALSE,group.intercept2=FALSE),
#'                                                    polynomial.degree=c(3,3))))
#'                                                    
#' Delta_AUC_2 <- Group_specific_Delta_AUC_estimation(MEM_Pol_group=MEM_estimation_2,
#'                                                    Group1="Group1",Group2="Group2",
#'                                                    time.G1=time_group1,time.G2=time_group2)
#'}

#' @seealso 
#'  \code{\link[AUCcomparison]{MEM_Polynomial_Group_structure}},
#'  \code{\link[AUCcomparison]{Group_specific_AUC_estimation}}
#' @rdname Group_specific_Delta_AUC_estimation
#' @export 
#' @import mathjaxr

Group_specific_Delta_AUC_estimation <- function(MEM_Pol_group,Group1,Group2,time.G1,time.G2,common.interval=TRUE,method="trapezoid",Averaged=FALSE){
  
  # Check of common.interval argument
  Check_commonTime <- ArgumentCheck::newArgCheck()
  if(isFALSE(is.logical(common.interval))){
    ArgumentCheck::addError(
      msg = "Error - The variable 'common.interval' must be a boolean",
      argcheck = Check_commonTime
    )
  }
  ArgumentCheck::finishArgCheck(Check_commonTime)
  
  
  if(common.interval == TRUE){
    min.time.interval <- max(min(time.G1,na.rm=TRUE),min(time.G2,na.rm=TRUE),na.rm=TRUE)
    max.time.interval <- min(max(time.G1,na.rm=TRUE),max(time.G2,na.rm=TRUE),na.rm=TRUE)
    
    if(min.time.interval > max.time.interval){
      stop("Impossible to estimate the difference of AUC on the common interval of time for the two groups: the inrestection of the two intervals is equal to an empty interval.")
    }else if(min.time.interval == max.time.interval){
      stop(paste("Impossible to estimate the difference of AUC on the common interval of time for the two groups: the inrestection of the two intervals is equal to a single time point:",min.time.interval))
    }else{
      time.G1 <- time.G1[which(time.G1 >= min.time.interval & time.G1 <= max.time.interval)]
      time.G2 <- time.G2[which(time.G2 >= min.time.interval & time.G2 <= max.time.interval)]
    }
  }
  
  AUC_Group1 <- Group_specific_AUC_estimation(MEM_Pol_group=MEM_Pol_group,time=time.G1,Groups=Group1,method=method,Averaged=Averaged) 
  AUC_Group2 <- Group_specific_AUC_estimation(MEM_Pol_group=MEM_Pol_group,time=time.G2,Groups=Group2,method=method,Averaged=Averaged) 
  Delta_AUC <- as.numeric(AUC_Group2 - AUC_Group1)
  return(Delta_AUC)
}