quasi_bayesian.R

# -----------------------------------------------------------
# Quasi-Bayesian estimation.
# Based on Słonecka et al. (2018), Radiation & Environ. Biophys. 57:195–203.
# -----------------------------------------------------------
#'
#' @param num_cases Number of different estimations to perform.
#' @param dics Number of observed dicentrics.
#' @param cells Number of total observed cells. 
#' @param coef_gamma Coefficients (C, β, γ) for gamma curve and std error.
#' @param coef_neutron Coefficients (C, α, 0) or (C, α, δ) for neutron curve and std error.
#' @param cov_gamma Covariance matrix for the Gamma curve. By default NULL.
#' @param cov_neutron Covariance matrix for the Neutron curve. By default NULL.
#' @param ratio_mu Approximate value (mean) for ratio.
#' @param ratio_sigma SD for ratio.
#' @param dist distribution for ratio. Can be "normal" or "gamma". By default "normal".
#' @param limr limits for the truncated normal distribution. e.g: c(0,3). By default NULL.
#' @example examples/qbayesian_dose_criticality_example.R
#' 
#' @import truncnorm numDeriv
#'
#' @return List containing estimated dose for gamma, neutron, 
#' and total, each with corresponding lower and upper 95% credibility interval 
#' bounds:
#'[[1]]
#'[[1]]$gamma
#'est     lwr     upr 
#'2.29400 2.08176 2.50624 
#'[[1]]$neutron
#'est       lwr       upr 
#'0.1990000 0.1297006 0.2682994 
#'[[1]]$total
#'est      lwr      upr 
#'2.493000 2.253929 2.732071 

qbayesian_dose_criticality <- function(num_cases, dics, cells, coef_gamma,
                                       coef_neutron, cov_gamma=NULL, 
                                       cov_neutron=NULL,ratio_mu, ratio_sigma,
                                       dist = "normal",
                                       limr = NULL) {
  
  
  #Check if all the needed information is provided:
  inputs <- list(num_cases, dics, cells, coef_gamma,
                 coef_neutron, ratio_mu, ratio_sigma,
                 dist)
  
  if (any(sapply(inputs, is.null))) {
    stop("Please, provide all the necessary inputs. Only limr, cov_gamma and cov_neutron can be NULL.")
  }
  
  if(is.null(limr)){
    limr <- c(0, 2 * ratio_mu)
  }
  
  #Calculate covariance matrix if is not provided:
  if(is.null(cov_gamma)){
    cov_gamma <- diag(coef_gamma[, "std.error"]^2)
  }
  
  if(is.null(cov_neutron)){
    cov_neutron <- diag(coef_neutron[, "std.error"]^2)
  }
  
  #If neutron curve is lineal, remember to put 0s:
  if (nrow(cov_gamma) != 3 || nrow(cov_neutron) != 3 ||
      nrow(coef_neutron) != 3 || nrow(coef_gamma) != 3) {
    stop("The input matrices must be of dimension 3: Provide the calibration 
          curve parameters and their covariance matrices for a quadratic model, 
          including the intercept (C). If the intercept or any quadratic term is 
          absent, it should be specified as zero")
  }
  
  
  #Select the estimates column:
  coef_gamma <- coef_gamma[, 1]
  coef_neutron <- coef_neutron[, 1]
  
  
  output <- list()
  
  for(i in 1:num_cases){
    #calculate yield and yield SD for each case:
    yf <- dics[[i]] / cells[[i]]
    sd.Y <- sqrt(dics[[i]]) / cells[[i]]
    
    #build the covariance matrix
    cov <- matrix(0, nrow = 8, ncol =  8)
    cov[1:3, 1:3] <- as.matrix(cov_gamma)
    cov[4:6, 4:6] <- as.matrix(cov_neutron)
    cov[7, 7] <- sd.Y^2
    cov[8, 8] <- ratio_sigma^2
    
    params <- c(coef_gamma, coef_neutron, yf, ratio_mu)
    
    dose_estimation_Q <- function(params){
      x1 <- params[1]
      x2 <- params[2]
      x3 <- params[3]
      
      x4 <- params[4]
      x5 <- params[5]
      x6 <- params[6]
      
      yf <- params[7]
      ratio_mu  <- params[8]
      
      # coefficients
      a <- x1 + x4/2
      b <- x2
      c <- x3
      d <- x5
      e <- x6
      
      # posterior for Dg
      P_Dg <- function(Dg){
        if(Dg <= 0){
          return(0)
        }
        
        r_est <- (yf - a - b*Dg - c*Dg^2) / (d*Dg)
        
        if(!is.finite(r_est) || r_est <= 0){
          return(0)
        }
        
        if(dist == "normal"){
          return(dtruncnorm(r_est, limr[1], limr[2], mean = ratio_mu, sd = ratio_sigma))
          
        }else{
          alpha_r <- (ratio_mu^2)/(ratio_sigma^2) 
          beta_r  <- (ratio_mu)/(ratio_sigma^2) 
          return(dgamma(r_est, shape = alpha_r, rate = beta_r))
        }
      }
      
      #maximize posterior
      opt <- optimize(f = P_Dg, interval = c(1e-6, 6), maximum = TRUE)
      Dg <- opt$maximum
      
      Dn <- tryCatch(
        {uniroot(function(Dn){
          e*Dn^2 + d*Dn + (a + b*Dg + c*Dg^2 - yf)
          }, c(0, 6))$root},
        error = function(err) {
          stop("Failed to find root for Dn: ", err$message)
      })
      
      return(c(Dg=Dg, Dn=Dn, Dt=Dg + Dn))
    }
    
    est <- dose_estimation_Q(params)
    
    est_gamma <- est[1]
    est_neutron <- est[2]
    est_total <- est[3]
    
    #Calculate SE
    J <- jacobian(func = dose_estimation_Q, x = params)
    var_est <- J %*% cov[1:8, 1:8] %*% t(J)  #Var(G(x)) = J*Var(x)*J^T
    se <- sqrt(diag(var_est))
    
    #Dose with upper and lower limits:
    est_gamma <- c(est = est_gamma, lwr = est_gamma - 1.96*se[1], upr = est_gamma + 1.96*se[1])
    est_neutron <- c(est = est_neutron, lwr = est_neutron - 1.96*se[2], upr = est_neutron + 1.96*se[2])
    est_total <- c(est = est_total, lwr = est_total - 1.96*se[3], upr = est_total + 1.96*se[3])
    names(est_gamma) <- names(est_neutron) <- names(est_total) <- c("est","lwr","upr")
    
    #output configuration for biodosetools:
    output[[i]] <- list(gamma = est_gamma, neutron = est_neutron, total = est_total)
    
  }
  return(output)
  
}