R/optimization.R

In nkurzaw/TPP2D: Detection of ligand-protein interactions from 2D thermal profiles (DLPTP)

## Internal functions for optimization of fits

.trim_sum <- function(x, na.rm = TRUE)
  # Trimmed sum (sum up all but the highest values)
{
  sum(x[-which(x == max(x, na.rm = TRUE))],
      na.rm = TRUE)
}

.min_RSS_h0 <- function(data, par, len_temp) {
  # Optimization function for fitting an intercept model to a protein's 2D
  # thermal profile by minimizing the sum of squared errors
  temp_i <- log2_value <- NULL
  
  beta_0 <- par[seq_len(len_temp)]
  
  sum(with(data, (beta_0[temp_i] - log2_value) ^ 2))
}

.min_RSS_h0_trim <- function(data, par, len_temp) {
  # Optimization function for fitting an intercept model to a protein's 2D
  # thermal profile by minimizing the trimmed sum of squared errors
  temp_i <- log2_value <- NULL
  
  beta_0 <- par[seq_len(len_temp)]
  
  .trim_sum(with(data, (beta_0[temp_i] - log2_value) ^ 2))
}

.min_RSS_h1 <- function(data, par, len_temp) {
  # Optimization function for fitting an dose-response model to a
  # protein's 2D thermal profile by minimizing the sum of squared errors
  temp_i <- log2_value <- log_conc <- NULL
  
  zeta <- par[1]
  slope <- par[2]
  beta_max <- par[3]
  beta_0 <- par[4:(len_temp + 3)]
  alpha <- par[(4 + len_temp):(3 + len_temp * 2)]
  
  sum(with(data, (
    beta_0[temp_i] + (alpha[temp_i] * beta_max) /
      (1 + exp(-slope * (log_conc - zeta))) -
      log2_value
  ) ^ 2))
}

.min_RSS_h1_slope_pEC50 <- function(data, par, len_temp) {
  # Optimization function for fitting an dose-response model to a
  # protein's 2D thermal profile by minimizing the sum of squared errors
  temp_i <- log2_value <- log_conc <- temperature <- NULL
  
  zeta <- par[1]
  zeta_slope <- par[2]
  slope <- par[3]
  beta_max <- par[4]
  beta_0 <- par[5:(len_temp + 4)]
  alpha <- par[(5 + len_temp):(4 + len_temp * 2)]
  
  sum(with(data,
           (
             beta_0[temp_i] + (alpha[temp_i] * beta_max) /
               (1 + exp(-slope * (
                 log_conc - (zeta + zeta_slope * temperature)
               ))) -
               log2_value
           )^ 2))
}


.min_RSS_h1_trim <- function(data, par, len_temp) {
  # Optimization function for fitting an dose-response model to a
  # protein's 2D thermal profile by minimizing the trimmed sum of
  # squared errors
  temp_i <- log2_value <- log_conc <- NULL
  
  zeta <- par[1]
  slope <- par[2]
  beta_max <- par[3]
  beta_0 <- par[4:(len_temp + 3)]
  alpha <- par[(4 + len_temp):(3 + len_temp * 2)]
  
  .trim_sum(with(data, (
    beta_0[temp_i] + (alpha[temp_i] * beta_max) /
      (1 + exp(-slope * (log_conc - zeta))) -
      log2_value
  ) ^ 2))
}

.min_RSS_h0_gradient <- function(data, par, len_temp) {
  # Analytically solved gradient function for .min_RSS_h1
  temp_i <- log2_value <- NULL
  
  beta_0 <- par[seq_len(len_temp)]
  
  data <-
    full_join(data,
              expand.grid(
                log_conc = unique(data$log_conc),
                temp_i = seq_len(max(data$temp_i))
              ),
              by = c("log_conc", "temp_i"))
  
  outer_dev <-
    with(data, 2 * (beta_0[temp_i] - log2_value))
  
  d_beta_0 <- apply(matrix(outer_dev, nrow = 5, byrow = TRUE),
                    2, sum, na.rm = TRUE)
  
  return(d_beta_0)
  
}

.min_RSS_h0_gradient_trim <- function(data, par, len_temp) {
  # Analytically solved gradient function for .min_RSS_h1
  temp_i <- log2_value <- NULL
  
  beta_0 <- par[seq_len(len_temp)]
  
  data <-
    full_join(data,
              expand.grid(
                log_conc = unique(data$log_conc),
                temp_i = seq_len(max(data$temp_i))
              ),
              by = c("log_conc", "temp_i"))
  
  outer_dev <-
    with(data, 2 * (beta_0[temp_i] - log2_value))
  
  d_beta_0 <- apply(matrix(outer_dev, nrow = 5, byrow = TRUE),
                    2, .trim_sum, na.rm = TRUE)
  
  return(d_beta_0)
  
}

.min_RSS_h1_gradient <- function(data, par, len_temp) {
  # Analytically solved gradient function for .min_RSS_h1
  temp_i <- log2_value <- log_conc <- NULL
  
  zeta <- par[1]
  slope <- par[2]
  beta_max <- par[3]
  beta_0 <- par[4:(len_temp + 3)]
  alpha <- par[(4 + len_temp):(3 + len_temp * 2)]
  
  data <-
    full_join(data,
              expand.grid(
                log_conc = unique(data$log_conc),
                temp_i = seq_len(max(data$temp_i))
              ),
              by = c("log_conc", "temp_i"))
  outer_dev <-
    with(data, 2 * (beta_0[temp_i] + (alpha[temp_i] * beta_max) /
                      (1 + exp(
                        -slope * (log_conc - zeta)
                      )) -
                      log2_value))
  d_zeta <- sum(outer_dev * with(data,-(alpha * beta_max) / (1 + exp(
    -slope * (log_conc - zeta)
  )) ^ 2 *
    exp(-slope * (log_conc - zeta)) * slope), na.rm = TRUE)
  d_slope <- sum(outer_dev * with(data,-(alpha * beta_max) / (1 + exp(
    -slope * (log_conc - zeta)
  )) ^ 2 *
    exp(-slope * (log_conc - zeta)) * (-(log_conc - zeta))), na.rm = TRUE)
  d_beta_max <- sum(outer_dev * with(data, alpha / (1 + exp(
    -slope * (log_conc - zeta)
  ))), na.rm = TRUE)
  d_beta_0 <-
    apply(matrix(outer_dev, nrow = 5, byrow = TRUE), 2, sum, na.rm = TRUE)
  d_alpha <- apply(matrix(
    outer_dev * with(data, beta_max / (1 + exp(
      -slope * (log_conc - zeta)
    ))),
    nrow = 5,
    byrow = TRUE
  ), 2, sum, na.rm = TRUE)
  return(c(d_zeta, d_slope, d_beta_max, d_beta_0, d_alpha))
  
}

.min_RSS_h1_gradient_trim <- function(data, par, len_temp) {
  # Analytically solved gradient function for .min_RSS_h1_trim
  temp_i <- log2_value <- log_conc <- NULL
  
  zeta <- par[1]
  slope <- par[2]
  beta_max <- par[3]
  beta_0 <- par[4:(len_temp + 3)]
  alpha <- par[(4 + len_temp):(3 + len_temp * 2)]
  
  data <-
    full_join(data,
              expand.grid(
                log_conc = unique(data$log_conc),
                temp_i = seq_len(max(data$temp_i))
              ),
              by = c("log_conc", "temp_i"))
  outer_dev <-
    with(data, 2 * (beta_0[temp_i] + (alpha[temp_i] * beta_max) /
                      (1 + exp(
                        -slope * (log_conc - zeta)
                      )) -
                      log2_value))
  d_zeta <- .trim_sum(outer_dev * with(data,-(alpha * beta_max) / (1 + exp(
    -slope * (log_conc - zeta)
  )) ^ 2 *
    exp(-slope * (log_conc - zeta)) * slope), na.rm = TRUE)
  d_slope <- .trim_sum(outer_dev * with(data,-(alpha * beta_max) / (1 + exp(
    -slope * (log_conc - zeta)
  )) ^ 2 *
    exp(-slope * (log_conc - zeta)) * (-(log_conc - zeta))), na.rm = TRUE)
  d_beta_max <- .trim_sum(outer_dev * with(data, alpha / (1 + exp(
    -slope * (log_conc - zeta)
  ))), na.rm = TRUE)
  d_beta_0 <-
    apply(matrix(outer_dev, nrow = 5, byrow = TRUE), 2, .trim_sum, na.rm = TRUE)
  d_alpha <- apply(matrix(
    outer_dev * with(data, beta_max / (1 + exp(
      -slope * (log_conc - zeta)
    ))),
    nrow = 5,
    byrow = TRUE
  ),
  2,
  .trim_sum,
  na.rm = TRUE)
  return(c(d_zeta, d_slope, d_beta_max, d_beta_0, d_alpha))
  
}

.min_RSS_h1_slope50_gradient <- function(data, par, len_temp) {
  # Analytically solved gradient function for .min_RSS_h1_slope_pEC50
  temp_i <- log2_value <- log_conc <- temperature <- NULL
  
  zeta <- par[1]
  zeta_slope <- par[2]
  slope <- par[3]
  beta_max <- par[4]
  beta_0 <- par[5:(len_temp + 4)]
  alpha <- par[(5 + len_temp):(4 + len_temp * 2)]
  
  data <-
    full_join(data,
              expand.grid(
                log_conc = unique(data$log_conc),
                temp_i = seq_len(max(data$temp_i))
              ),
              by = c("log_conc", "temp_i"))
  outer_dev <-
    with(data,
         2 * (beta_0[temp_i] + (alpha[temp_i] * beta_max) /
                (1 + exp(
                  -slope * (log_conc - (zeta + zeta_slope * temperature))
                )) -
                log2_value))
  
  d_zeta <- sum(outer_dev *
                  with(data,-(alpha * beta_max) /
                         (1 + exp(
                           -slope * (log_conc - (zeta + zeta_slope * temperature))
                         )) ^ 2 *
                         exp(-slope * (
                           log_conc - (zeta + zeta_slope * temperature)
                         )) * slope),
                na.rm = TRUE)
  
  d_zeta_slope <- sum(outer_dev * with(
    data,
    -(alpha * beta_max) / (1 + exp(-slope * (
      log_conc - (zeta + zeta_slope * temperature)
    ))) ^ 2 *
      exp(-slope * (
        log_conc - (zeta + zeta_slope * temperature)
      )) * (slope * temperature)
  ), na.rm = TRUE)
  
  d_slope <- sum(outer_dev * with(data,-(alpha * beta_max) / (1 + exp(
    -slope * (log_conc - (zeta + zeta_slope * temperature))
  )) ^ 2 *
    exp(-slope * (
      log_conc - (zeta + zeta_slope * temperature)
    )) * (-(
      log_conc - (zeta + zeta_slope * temperature)
    ))), na.rm = TRUE)
  d_beta_max <- sum(outer_dev * with(data, alpha / (1 + exp(
    -slope * (log_conc - (zeta + zeta_slope * temperature))
  ))), na.rm = TRUE)
  d_beta_0 <-
    apply(matrix(outer_dev, nrow = 5, byrow = TRUE), 2, sum, na.rm = TRUE)
  d_alpha <- apply(matrix(
    outer_dev * with(data, beta_max / (1 + exp(
      -slope * (log_conc - (zeta + zeta_slope * temperature))
    ))),
    nrow = 5,
    byrow = TRUE
  ), 2, sum, na.rm = TRUE)
  return(c(d_zeta, d_zeta_slope, d_slope, d_beta_max, d_beta_0, d_alpha))
  
}