R/ccvx_fit.R

In paulmanser/concavex: Concavex Bayesian Dose-Response Modeling

#' Sample from Concavex model posteriors
#' 
#' This function performs Gibbs sampling from posterior densities of concavex model parameters as well as transformations of parameters.
#' 
#' @param ccvx.mod  JAGS model file as specified by \link{ccvx_build_jags}
#' @param doses A vector of dose strengths, with placebo listed first
#' @param mu.hat A vector of parameter estimates for each of the doses
#' @param std.err A vector of standard errors for each of the doses
#' @param n.chains Number of chains used to for Gibbs sampling. Default is 4
#' @param gibbs.samples Number of samples to draw for each chain after burn in. Default is 5000
#' @param burn.in Number of burn in samples to be drawn. Default is 1000
#' @param sd.ph3 The standard deviation of the endpoint being assessed in phase 3. Only needs to be specified when ccvx_fit() computes posterior predictive probabilities (i.e. model-based DDCPs)
#' @param n.per.arm.ph3 The number of patients per arm being assessed in hypothetical phase 3 study. Only needs to be specified when ccvx_fit() computes posterior predictive probabilities (i.e. model-based DDCPs)
#' 
#' @return A list with the elements
#' \item{ccvx.mod}{The JAGS code used for Gibbs sampling as generated by \link{ccvx_build_jags} or \link{ccvx5_build_jags}}
#' \item{jags.samples}{output containing Gibbs samples from parameter posteriors generated}
#' \item{coda.samples}{output containing coda samples for MCMC diagnostics}
#' \item{doses}{A vector of dose strengths, with placebo listed first}
#' \item{mu.hat}{A vector of parameter estimates for each of the doses}
#' \item{std.err}{A vector of standard errors for each of the doses}
#' 
#' @export
#' @examples 
#' ccvx.mod <- ccvx_build_jags()
#' ccvx.samples <- ccvx_fit(ccvx.mod, doses = 0:4, mu.hat = c(1, 20, 50, 60, 65), std.err = rep(10, 5))
#' names(ccvx.samples$jags.samples)
#' ccvx_plot_fit(ccvx.samples)

ccvx_fit <- function(ccvx.mod, doses, mu.hat, std.err, n.chains = 4, gibbs.samples = 5000, burn.in = 1000,
                     sd.ph3 = NULL, n.per.arm.ph3 = NULL) {
  
  if(any(grep('trt.post.pred.dose', ccvx.mod)) & (is.null(sd.ph3) | is.null(n.per.arm.ph3)  ))
    stop("The 'ccvx.mod' JAGS script has been specified to compute posterior predictive probabilities. \n Please provide values for 'n.per.arm.ph3' and 'sd.ph3.'")
  
  if(any(grep('5-parameter concavex spline model', ccvx.mod))){ 
    post.vars <- c("theta_0", "theta_1", "gamma_1", "gamma_2", "alpha",
                      "mu.tilde", "dose.post", "trt.post", "trt.post.dose")
    
    post.pred.vars <- c("theta_0", "theta_1", "gamma_1", "gamma_2", "alpha", "mu.tilde", "dose.post", "trt.post", 
                           "trt.post.dose", "mu.tilde.pred", "trt.post.pred", "mu.tilde.pred.dose", "trt.post.pred.dose")
    
    coda.vars <- c("theta_0", "theta_1", "gamma_1", "gamma_2", "alpha")
  }
  if(any(grep('3-parameter concavex model', ccvx.mod))){ 
    post.vars <- c("theta_0", "theta_1", "lambda", "mu.tilde", "dose.post", "trt.post", "trt.post.dose", "ed90")
    
    post.pred.vars <- c("theta_0", "theta_1", "lambda", "mu.tilde", "dose.post", "trt.post", "trt.post.dose", "ed90",
                           "mu.tilde.pred", "trt.post.pred", "mu.tilde.pred.dose", "trt.post.pred.dose")
    
    coda.vars <- c("theta_0", "theta_1", "lambda")
  }
  
  tc.ccvx <- textConnection(ccvx.mod)

  if(!any(grep('trt.post.pred.dose', ccvx.mod))) {
  # initialize model
  model.init <- jags.model(tc.ccvx,
                           data = list('dose' = doses / max(doses),
                                       'eff' = mu.hat,
                                       'tau' = 1 / std.err^2,
                                       'pred.doses' = seq(0, 1, length.out = 250)),
                           n.chains = n.chains,
                           n.adapt = burn.in)
  # gibbs sampling
  ccvx.samples <- jags.samples(model.init, post.vars, n.iter = gibbs.samples)
  
  }else{
    std.err.ph3 <- sd.ph3 / sqrt(n.per.arm.ph3)
    # initialize model
    model.init <- jags.model(tc.ccvx,
                             data = list('dose' = doses / max(doses),
                                         'eff' = mu.hat,
                                         'tau' = 1 / std.err^2,
                                         'pred.doses' = seq(0, 1, length.out = 250),
                                         'tau.ph3' = 1 / std.err.ph3^2),
                             n.chains = n.chains,
                             n.adapt = burn.in)
    
    # gibbs sampling
    ccvx.samples <- jags.samples(model.init, post.pred.vars, n.iter = gibbs.samples)
  }
  
  close(tc.ccvx)
  
  # diagnostics with coda
  coda.samples <- coda.samples(model.init, coda.vars, gibbs.samples)
  
  out <- list(ccvx.mod = ccvx.mod, jags.samples = ccvx.samples, coda.samples = coda.samples,
              doses = doses, mu.hat = mu.hat, std.err = std.err)
  out
}