R/SimPatientOutcomes.MVNByGroup.R

In kwathen/OCTOPUS: OCTOPUS

##### COPYRIGHT #############################################################################################################
#
# Copyright (C) 2018 JANSSEN RESEARCH & DEVELOPMENT, LLC
# This package is governed by the JRD OCTOPUS License, which is the
# GNU General Public License V3 with additional terms. The precise license terms are located in the files
# LICENSE and GPL.
#
#############################################################################################################################.

#Formerley called SimPatientOutcomes.MVNWithCovariate but name was changed to be more descriptive
#####################################################################################
#' @name SimPatientOutcomes.MVNByGroup
#' @aliases{SimPatientOutcomes.MVNByGroup}
#' @title SimPatientOutcomes.MVNByGroup
#' @description  In this SimPatientOutcomes function we simulate patients from a MVN where there is a group covariate.
#'   For each arm in the ISA you specify a lSimArm that specifies how the patients in that
#'   arm are to be simulated.
#'   The cSimoutcomes is a structure( list(lSimArm1, lSimArm2,..., lSimArmXX), class="MVNByGroup" )
#'   Assuming the number of groups you want to simulate from is J, each lSimArm will have the following elements
#'
#'   Required for XX in 1:J :
#'       $vMeanXX - The mean vector for group XX
#'       $mVarCovXX - The var-cov matrix for group XX
#'
#'       $vProbGroup - Probability a patient is in each group
#'       $vObsTime   - The times each of the outcomes are observed length( vObsTime ) == length( vMeanXX ); for all XX
#'
#'   Optional Argument to constrain the simulated results if needed:
#'       $dMinimum - All simulated values must be >= dMinimum
#'       $dMaximum - All simulated values must be <= dMaximum
#'
#
#' @examples{ ## This example contains 2 arms and each arm will simulate from 2 populaitons.
#' ## Assuming vMeanR1, vMeanNR2 are mean vector for arm 1 defined elsewhere in the code and have length = 3
#' ## Assuming mVarCovR1, mVarCovNR1 are var-cov for arm 1 defined elsewhere in the code and are a 3x3 matrix
#' ## Assuming vMeanR2, vMeanNR2 are mean vector for arm 2 defined elsewhere in the code and have length = 3
#' ## Assuming mVarCovR2, mVarCovNR2 are var-cov for arm 2 defined elsewhere in the code and are a 3x3 matrix
#'
#' ## The following example would simulate patient with 3 measurements.  For arm 1, on average, 20% would
#' ## be from group 1 and 80% from group 2 with the respective mean and var-cov.  Patients in arm 2 would be simulated
#' ## with 30%, on average, from group 1 and 70%, on average, from group 2 with the respective mean and var cov.
#' vMeanR1    <- c( 10, 15 )
#' mVarCovR1  <- diag(1, nrow=2)
#' vMeanNR1   <- c( 5, 8 )
#' mVarCovNR1 <- diag( 1, nrow=2)
#' lSimArm1   <- list( vMean1     = vMeanR1,
#'                     mVarCov1   = mVarCovR1,
#'                     vMean2     = vMeanNR1,
#'                     mVarCov2   = mVarCovNR1,
#'                     vProbGroup = c( 0.2, .8 ),
#'                     dMinimum   = 0,
#'                     vObsTime   = c( 0, 2, 4 ))
#'
#' vMeanR2    <- c( 10, 15 )
#' mVarCovR2  <- diag(1, nrow=2)
#' vMeanNR2   <- c( 5, 8 )
#' mVarCovNR2 <- diag( 1, nrow=2)
#' lSimArm2   <- list( vMean1     = vMeanR2,
#'                     mVarCov1   = mVarCovR2,
#'                     vMean2     = vMeanNR2,
#'                     mVarCov2   = mVarCovNR2,
#'                     vProbGroup = c( 0.3, .7 ),
#'                     dMinimum   = 0,
#'                     vObsTime   = c( 0, 2, 4 ))
#'
#'  cSimOutcomes <- structure( list(lSimArm1 = lSimArm1, lSimArm2 = lSimArm2), class =  "MVNByGroup" )
#' }
#' @seealso { \href{https://github.com/kwathen/OCTOPUS/blob/master/R/SimPatientOutcomes.MVNByGroup.R}{View Code on GitHub} }
#' @export
SimPatientOutcomes.MVNByGroup <- function( cSimOutcomes,  cISADesign , dfPatCovISA )
{
    if( !is.null(  dfPatCovISA  ) )
        stop( "SimPatientOutcomes.MVNByGroup is not designed to incorporate patient covariates and  dfPatCovISA  is not NULL, please try SimPatOutcomes.MVNWithCovariates.")
    #print("SimPatientOutcomes.MVNByGroup" )

    vQtyPats        <- cISADesign$vQtyPats
    nQtyGroups      <- length( cSimOutcomes$lSimArm1$vProbGroup)
    vPatTrt         <- rep( cISADesign$vTrtLab, vQtyPats )
    nQtyArms        <- length( vQtyPats )

    vstrMeanName    <- paste( "vMean", 1:nQtyGroups, sep="")
    vstrVarCovName  <- paste( "mVarCov", 1:nQtyGroups, sep="")
    vstrSimArm      <- paste( "lSimArm", 1:nQtyArms, sep="" )

    ###################################################################
    #TODO: Validation to be moved to a Validate phase?
    bValid <- TRUE
    strErr <- ""
    strErrPrfx <- "SimPatientOutcomes.MVNByGroup - Validation Error: "
    if( length( cSimOutcomes ) != nQtyArms )
    {
        bValid <- FALSE
        strErr <- paste( strErrPrfx, "The number of arms in ISA ( ", nQtyArms, " )")
        strErr <- paste( strErr, " != number of lSimArms ( ", length( cSimOutcomes), " in cSimOutcomes object.")
    }

    for( iArm in 1:nQtyArms )
    {
        lSimArm <- cSimOutcomes[[ vstrSimArm[ iArm ]]]
        if( sum( lSimArm$vProbGroup )  != 1 )
        {
            bValid <- FALSE
            strErr <- paste( strErr, strErrPrfx,  "For lSimArm", iArm, " sum( vProbGroup ) != 1.")

        }
        if( length( lSimArm$vMean1 ) != length( lSimArm$vObsTime ) )
        {
            bValid <- FALSE
            strErr <- paste( strErr, strErrPrfx,  "For lSimArm", iArm, "  length( lSimArm$vMean1 ) != legnth( lSimArm$vObsTime )")
        }

        #TODO: More validation
        # 1.  Make sure the number of vMean and mVarCov supplied = length( vProbGroup )
        # 2.  Make sure length( vMeanX ) = nrow( mVarCovX ) = ncol( mVarCov )
    }

    ### End Validation
    ###################################################################

    mOutcome        <- NULL
    vIndGrp         <- NULL   #Will be a vector of what group a patient is in


    iGrp <- 1
    iArm <- 1
    repeat  # loop over the arms
    {

        nQtyPats    <- vQtyPats[ iArm ]    # This is the number of patients on the arm, now need to know how many are IR and NV
        lSimArm     <- cSimOutcomes[[ vstrSimArm[ iArm ]]]

        vProbGroup  <- lSimArm$vProbGroup
        vPatsPerGrp <- rmultinom( 1, size= nQtyPats, prob = vProbGroup )[,1]  #rmultinom returns a matrix with column for each sample, only need 1 column
        vIndGrp     <- c( vIndGrp, rep( 1:nQtyGroups, vPatsPerGrp) )

        # Now need to simulate outcomes from each group
        iGrp       <- 1

        repeat  #Loop over the groups
        {
            vMean      <- lSimArm[[ vstrMeanName[ iGrp ] ]]
            mVarCov    <- lSimArm[[ vstrVarCovName[ iGrp ] ]]  #var-cov only depends on the group

            if(vPatsPerGrp[ iGrp ] > 0 )
            {
                mOutcomeTmp <-  MASS::mvrnorm( n= vPatsPerGrp[ iGrp ], vMean, mVarCov, tol = .1 )

                ## Need to check the minimums and maximums
                if( is.null( lSimArm$dMinimum ) == FALSE )  # A minimum was supplied
                {
                    mOutcomeTmp[ mOutcomeTmp <= lSimArm$dMinimum  ] <- lSimArm$dMinimum
                }

                if( is.null( lSimArm$dMaximum ) == FALSE )  # A maximum was supplied
                {
                    mOutcomeTmp[ mOutcomeTmp >= lSimArm$dMaximum  ] <- lSimArm$dMaximum
                }
                mOutcome    <- rbind( mOutcome, mOutcomeTmp)
            }

            if( iGrp == nQtyGroups )
                break
            iGrp <- iGrp + 1
        } #End of for loop over the groups

        if( iArm == nQtyArms )
            break

        iArm <- iArm + 1

    } #End of loop over the arms

    lSimDataRet <- structure( list( mSimOut1 = mOutcome, vObsTime1 = cSimOutcomes[[1]]$vObsTime, vCov1 = vIndGrp  ), class= class(cSimOutcomes) )

    lSimDataRet$nQtyOut  <- 1
    lSimDataRet$vPatTrt  <- vPatTrt
    #lSimDataRet$vPatISA  <- rep( cTrialDesign$vISA, cTrialDesign$vQtyPats )
    return( lSimDataRet )
}