analytical.R

# -----------------------------------------------------------
# Analytical estimation.
# Based on Słonecka et al. (2018), Radiation & Environ. Biophys. 57:195–203.
# Based on an original script from Biodose Tools by David Endesfelder.
# -----------------------------------------------------------
#'
#' @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. If NULL, an approximation will be calculated.
#' @param cov_neutron Covariance matrix for the Neutron curve. By default NULL. If NULL, an approximation will be calculated.
#' @param ratio neutron/gamma dose ratio ρ.
#' @example examples/analytical_dose_criticality_example.R
#' 
#' @import numDeriv
#'
#' @return List containing estimated dose for gamma, neutron, 
#' and total, each with corresponding lower and upper 95% confidence 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 

analytical_dose_criticality <- function(num_cases, dics, cells, coef_gamma, 
                                        coef_neutron, cov_gamma = NULL,
                                        cov_neutron = NULL, ratio){
                                       
  #Check if all the needed information is provided:
  inputs <- list(num_cases, dics, cells, coef_gamma,
                 coef_neutron, ratio)
  
  if (any(sapply(inputs, is.null))) {
    stop("Please, provide all the necessary inputs. Only cov_gamma and cov_neutron can be NULL.")
  }
  
  #Calculate covariance matrix if is not provided - just approximation:
  if(is.null(cov_gamma)){
    cov_gamma <- diag(coef_gamma[, "std.error"]^2)
    rownames(cov_gamma) <- colnames(cov_gamma) <- rownames(coef_gamma)
  }
  
  if(is.null(cov_neutron)){
    cov_neutron <- diag(coef_neutron[, "std.error"]^2)
    rownames(cov_neutron) <- colnames(cov_neutron) <- rownames(coef_neutron)
  }
  
  #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()
  
  #iterate for each case:
  for(i in 1:num_cases){
    r <- ratio
    #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 = 7, ncol =  7)
    cov[1:3, 1:3] <- as.matrix(cov_gamma)
    cov[4:6, 4:6] <- as.matrix(cov_neutron)
    cov[7, 7] <- sd.Y^2
    
    #coefficients: yf = C + alpha*Dn + beta*Dg + gamma*Dg^2 + delta*Dn^2
    # C = (coef_neutron[1] + coef_gamma[1])/2
    # alpha = coef_neutron[2]
    # beta = coef_gamma[2]
    # gamma = coef_gamma[3]
    # delta = coef_neutron[3] = tmp
    # yf = dics/cells
    
    # a*Dg^2 + b*Dg + c = 0 :
    tmp <- coef_neutron[3] 
    a <- coef_gamma[3] + ((ifelse(tmp == 0, 0, tmp))*r^2)
    b <- coef_gamma[2] + (coef_neutron[2]*r)
    c <- ((coef_neutron[1] + coef_gamma[1])/2) - yf
    
    #dose estimations:
    est_gamma <- (-b + sqrt(b^2-4*a*c))/(2*a)
    est_neutron <- r * est_gamma
    est_total <- est_gamma + est_neutron
    
    #x1 = coef_gamma[1]
    #x2 = coef_gamma[2] 
    #x3 = coef_gamma[3] 
    #x4 = coef_neutron[1]
    #x5 = coef_neutron[2] 
    #x6 = coef_neutron[3] = 0 (linear) != 0 (quadratic)
    #x7 = yf 
    
    #Dg = (-b + sqrt(b^2-4*a*c))/(2*a); 
    
    #define the parameters for the dose estimation function:
    params <- c(coef_gamma, coef_neutron, yf) 
    
    #dose estimation using analytical method:
    dose_estimation_A <- function(params) {#function G(X) --> Var(G(x)) = J*Var(x)*J^T
      x1 <- params[1]; x2 <- params[2]; x3 <- params[3]; 
      x4 <- params[4]; x5 <- params[5]; x6 <- params[6]; x7 <- params[7]
      
      if(tmp != 0){ # quadratic neutron curve
        Dg <- (-(x2 + x5*r) + sqrt((x2 + x5*r)^2 - 4*(x3 + x6*r^2)*((x4 + x1)/2 - x7))) / (2*(x3 + x6*r^2))
      } else {      # linear neutron curve
        Dg <- (-(x2 + x5*r) + sqrt((x2 + x5*r)^2 - 4*x3*((x4 + x1)/2 - x7))) / (2*x3)
      }
      
      Dn <- Dg * r
      
      return(c(Dg=Dg, Dn=Dn, Dtot=Dg+Dn))
    }
    
    #Calculate SE
    J <- jacobian(func = dose_estimation_A, x = params)
    var_est <- J %*% cov[1:7, 1:7] %*% 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)
  
}