R/optimization.R

In nkurzaw/TPP2Drep: FDR-controlled analysis of 2D-TPP experiments with replicates

## 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_rep <- function(data, par, len_rep_temp) {
    # Optimization function for fitting an intercept model to a protein's 2D
    # thermal profile by minimizing the sum of squared errors
    rep_temp_i <- log2_value <- NULL
    #len_rep_temp <- max(data$rep_temp_i)
    beta_0 <- par[seq_len(len_rep_temp)]
    
    sum(with(data, (beta_0[rep_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_slope_pEC50_rep <- function(data, par, len_temp, len_rep_temp){
    # Optimization function for fitting an dose-response model to 
    # replicates of 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_rep_temp + 4)]
    alpha <- par[(len_rep_temp + 5):(len_rep_temp + len_temp + 4)]
    
    with(data,
         sum((beta_0[rep_temp_i] + (alpha[temp_i] * beta_max)/
                  (1 + exp(-slope * 
                               (log_conc - 
                                    (zeta + zeta_slope * temperature)))) - 
                  log2_value)^2))
    
}
# min_RSS_h1_cpp <- 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
#   # using a CPP implementation of the optimization task
# {
#   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, rcpp_compute_residuals(temp_i, zeta, slope, beta_max,
#                                       beta_0, alpha, log_conc, log2_value))
#   )
# }

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_h1_trim_cpp <- 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 using a CPP implementation of the optimization
#   # task
# {
#   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, rcpp_compute_residuals(temp_i = temp_i, zeta = zeta,
#                                       slope = slope, beta_max = beta_max,
#                                       beta_0 = beta_0, alpha = alpha,
#                                       log_conc = log_conc,
#                                       log2_value = log2_value))
#   )
# }

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))
  
}

min_RSS_h1_slope50_rep_gradient <- function(data, par, len_temp) {
    # Analytically solved gradient function for min_RSS_h1_slope_pEC50
    temp_i <- log2_value <- log_conc <- temperature <- NULL
    unique_rep <- unique(data$replicate)
    len_rep <- length(unique_rep)
    len_rep_temp <- max(data$rep_temp_i)
    unique_temp <- unique(data$temperature)
    
    zeta <- par[1]
    zeta_slope <- par[2]
    slope <- par[3]
    beta_max <- par[4]
    beta_0 <- par[5:(len_rep_temp + 4)]
    alpha <- par[(len_rep_temp + 5):(len_rep_temp + len_temp + 4)]
    
    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[rep_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[temp_i] * 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[temp_i] * 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[temp_i] * 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[temp_i] / (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))
    
}