R/eff_summary.R

Defines functions eff_summary eff_suggest

Documented in eff_suggest eff_summary

eff_summary <- function(eff.structure, eff.Sigma, eff.sd_trans, n.sim = 300000,
                        seed = 123, plot.flag = F, plot.title = T){

  #' Compute the summary statistics of efficacy measure with specified
  #' parameters.
  #'
  #' Numerically compute the Mean, mariginal standard deviance and mariginal
  #' correlation matrix of efficacy measure generated with specified parameters.
  #' The function provides plots of marginal density of generated efficacy and
  #' correlation matrix for each dose. Check details to see the efficacy data
  #' generation procedures.
  #'
  #' @param eff.structure  A matrix providing the mean of the multivariate
  #'   Gaussian distribution in efficacy data generation. Specifically, the
  #'   \eqn{(i, j)}th element represents the mean value of \eqn{i}th dose and
  #'   \eqn{j}th cycle of the Gaussian distribution for efficacy data
  #'   generation.
  #' @param eff.Sigma The covariance matrix of the multivariate Guassian
  #'   distribution in efficacy data generation.
  #' @param eff.sd_trans A positive number controling the skewness of the
  #'   distribution of the efficacy response.
  #' @param n.sim Number of simulations for calculation summary statistics.
  #'   Default is 300,000
  #' @param seed 	The seed of R's random number generator. Default is 123
  #' @param plot.flag Whether output the marginal density, and correlation
  #'   matrix or not. Default is FALSE.
  #' @param plot.title Whether display the title of the plot or not. Default is
  #'   TRUE
  #'
  #' @details The user can simulate longitudinal efficacy response with
  #'   different dose-efficacy and cycle-efficacy pattern using argument
  #'   \code{eff.structure}, \code{eff.Sigma} and \code{eff.sd_trans}. The
  #'   sampling process of efficacy response starts from generating \eqn{z =
  #'   {z1, \ldots, zd} } from multivariate Gaussian distribution \deqn{z ~
  #'   MVN(\mu, V)}, where \eqn{\mu} and \eqn{V} are specified by
  #'   \code{eff.structure} and \code{eff.Sigma}, respectively. Define
  #'   \eqn{\phi} be the density of \eqn{N(0, \sigma^2)} with CDF \eqn{\Phi},
  #'   where \eqn{\sigma^2} is set by \code{eff.sd_trans}. Then the efficacy
  #'   measure is generated by taking the CDF of \eqn{z}: \deqn{x={x1, \ldots,
  #'   xd} = \Phi(z) = { \Phi(z1), \ldots, \Phi(zd)}}. Notice here the variance
  #'   parameter \eqn{\sigma^2_{trans}} controls the variance of the generated
  #'   efficacy.
  #'
  #' @return \item{eff.M}{A matrix recording the efficacy mean whose \eqn{(i,
  #' j)}th element represents the efficacy mean of \eqn{i}th dose level and
  #' \eqn{j}th cycle }
  #'
  #' \item{eff.cor.ls}{A list with a length of dose levels numbers recording the
  #' marginal correlation matrix across cycles of efficacy data for each dose
  #' level}
  #'
  #' @examples
  #' data("eff")       # load eff.RData from package phase1PRMD. Details see "?eff"
  #' eff.structure = eff$Dose_Cycle_Meff["plat", "dec", , ]
  #' eff.Sigma = eff$Sigma
  #' eff.sd_trans = eff$sd_trans
  #'
  #'# res <- eff_summary(eff.structure, eff.Sigma, eff.sd_trans, n.sim = 300000,
  #'#                    seed = 123)
  #'# res
  #' # set a special cases and check the density and correlation plots
  #'# eff_summary(eff.structure = matrix(eff.structure[cbind(c(1:6), c(1:6))],
  #'#                                   nrow = 1, ncol = 6),
  #'#             eff.Sigma, eff.sd_trans, n.sim = 300000, seed = 123,
  #'#             plot.flag = TRUE, plot.title = FALSE)
  #'
  #' @import ggplot2
  #' @importFrom gridExtra grid.arrange
  #' @importFrom MASS mvrnorm
  #' @importFrom reshape2 melt
  #' @export

  set.seed(seed)
  n.doses <- dim(eff.structure)[1]
  n.cycles <- dim(eff.structure)[2]
  eff.M <- matrix(NA, nrow = n.doses, ncol = n.cycles)
  rownames(eff.M) <- paste0("D", 1:n.doses)
  colnames(eff.M) <- paste0("c", 1:n.cycles)
  cormat.ls <- list()
  eff.cor.ls <- list()
  dose_cycle_eff_d <- NULL
  for(d in 1:n.doses){
    z <- mvrnorm(n.sim, mu = eff.structure[d, ], Sigma = eff.Sigma)
    # the efficacy data is the cdf of the correlated data from nultivariate normal
    eff.gen <- pnorm(z, mean = rep(0, n.doses), sd = eff.sd_trans)
    eff.gen.sd <- format(round(sqrt(diag(var(eff.gen))), 2), 2)
    colnames(eff.gen) <- paste0("c", 1:n.cycles, " sd:", eff.gen.sd)
    eff.M[d, ] <- colMeans(eff.gen)
    cormat <- round(cor(eff.gen), 2)
    eff.cor.ls[[d]] <- cormat
    if(plot.flag == T){
      cormat[upper.tri(cormat)] <- NA
      cormat.ls[[d]] <- melt(cormat, na.rm = TRUE)
      cormat.ls[[d]]$text <- cormat.ls[[d]]$value
      cormat.ls[[d]]$text[which(cormat.ls[[d]]$Var1 == cormat.ls[[d]]$Var2)] <-
        paste0("sd:", eff.gen.sd)
      dose_cycle_eff_d <- data.frame(eff = as.vector(t(eff.gen)),
                                     cycle = rep(1:n.cycles, n.sim))
      dose_cycle_eff_d$cycle <- factor(dose_cycle_eff_d$cycle,
                                       levels = c(1:n.cycles),
                                       labels = paste("cycle", 1:n.cycles))
      # plot the efficacy density and correlation matrix
      pd <- ggplot() +
        geom_density(aes(x = dose_cycle_eff_d$eff, color = cycle,
                         linetype = cycle, fill = cycle),
                     data = dose_cycle_eff_d, alpha = 0.3) + theme_bw() +
        theme(axis.title.x = element_blank(), axis.title.y = element_blank())
      pc <- ggplot(data = cormat.ls[[d]], aes(cormat.ls[[d]]$Var1,
                                              cormat.ls[[d]]$Var2,
                                              fill = cormat.ls[[d]]$value))+
        geom_tile(color = "white") +
        scale_fill_gradient2(low = "blue", high = "purple", mid = "white",
                             midpoint = 0, limit = c(-1, 1), space = "Lab",
                             name = element_blank()) +
        theme_minimal()+ # minimal theme
        coord_fixed() +
        geom_text(aes(cormat.ls[[d]]$Var1, cormat.ls[[d]]$Var2,
                      label = cormat.ls[[d]]$value), color = "black",
                  size = 3) +
        theme(
          axis.title.x = element_blank(),
          axis.text.x = element_blank(),
          axis.title.y = element_blank(),
          panel.grid.major = element_blank(),
          panel.border = element_blank(),
          panel.background = element_blank(),
          axis.ticks = element_blank()) +
        guides(fill=FALSE)
      print(grid.arrange(pd, pc, widths = c(4.5, 2.5), ncol = 2,
                         top = ifelse(plot.title == T,
                                      paste("Efficacy data summary, Dose ", d),
                                      "")))
      readline(prompt="Press [enter] to check the next plot")
    }
  }
  return(list(eff.M = eff.M, eff.cor.ls = eff.cor.ls))
}



eff_suggest <- function(eff.M, eff.sd, eff.sd_trans, n.sim = 30000){
  #' Suggest the input \code{eff.structure} of function \code{SimPRMD} with
  #' selected \code{eff.sd_tran}
  #'
  #' Suggest the input \code{eff.structure} of function SimPRMD  with selected
  #' \code{eff.sd_tran} for given efficacy mean matrix and efficacy standard
  #' deviation
  #'
  #' @param eff.M    The efficacy mean matrix whose (i,j)th element save the
  #'   target mean of the efficacy data
  #' @param eff.sd   The target standard deviation matrix for all dose and
  #'   cycles. Notice that the generated efficacy may have different standard
  #'   deviation matrix due to the correlations across cycles
  #' @param eff.sd_trans The eff.sd_trans for test. Notice variance of the
  #'   generated efficacy data will be effected by \code{eff.sd_trans}.
  #' @param n.sim    The number of simulations for the numerical calculation in
  #'   the function. The default is 30,000
  #'
  #' @return \item{eff.suggest}{The matrix suggested for the input
  #'   \code{eff.structure} of function \code{SimPRMD}}
  #'
  #' @importFrom stats rbeta qnorm
  #' @export
  #'
  #'
  #' @examples
  #' # Provide an target efficacy mean matrix for all dose and cycles
  #' eff.M <- matrix(rep(3:8/10, 6), nrow = 6, ncol = 6)
  #'
  #' # Give a target standard deviation matrix for all dose and cycles
  #' # Notice that the generated efficacy may have difference standard deviation
  #' # matrix due to the correlations across cycles
  #' eff.sd <- matrix(0.2, nrow = 6, ncol = 6)
  #'
  #' # Select a eff.sd_trans for testing. The efficacy variance are mainly
  #' # controlled by the eff.sd_trans
  #' eff.sd_trans <- 1.5  # or other positive value
  #' eff.structure <- eff_suggest(eff.M = eff.M, eff.sd = eff.sd,
  #'                              eff.sd_trans = eff.sd_trans)
  #'
  #' # check whether the suggested eff.M and the selected sd_trans
  #' # generate the desirable scenario
  #' eff.Sigma <- diag(6)
  #' diag(eff.Sigma[-1,]) = 0.5
  #' diag(eff.Sigma[, -1]) = 0.5
  #' eff.check <- eff_summary(eff.structure = eff.structure,
  #'                          eff.Sigma = eff.Sigma,
  #'                          eff.sd_trans = eff.sd_trans,
  #'                          plot.flag = FALSE)
  #' eff.check$eff.M
  #' eff.check$eff.cor.ls
  #'
  #'

  n.doses <- nrow(eff.M)
  n.cycles <- ncol(eff.M)
  eff.suggest <-  sapply(1:n.cycles, function(c){
    sapply(1:n.doses, function(d){
      m = eff.M[d, c]
      v <- eff.sd[d, c]^2
      sim.eff <- rbeta(n.sim, shape1 = m*(m*(1-m)/v - 1),
                       shape2 = (1-m)*(m*(1-m)/v - 1))
      sim.eff.z <- qnorm(sim.eff, mean = 0, sd = eff.sd_trans)
      mean(sim.eff.z)
    })
  })
  eff.suggest
}
LuZhangstat/Phase1Mayo documentation built on March 16, 2020, 1:32 a.m.